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PREFACE 


This  report  describes  methods  for  the  design  and  Implementation  of 
digital  filters  that  have  application  in  modern  aircraft  avionics  and 
flight  control  systems.  The  first  section  of  the  report  deals  with  the 
design  of  digital  filters  and  the  interaction  between  the  design  and 
implementation  tasks.  The  second  section  describes  the  implementation 
of  digital  filters  and  the  procedures  that  are  used  to  determine  the 
computer  word  length.  In  the  third  section  the  use  of  multiplexing 
in  the  implementation  of  digital  filters  is  described.  The  last  section 
is  concerned  with  conclusions  and  recommendations  for  future  work  that 
will  result  in  less  expensive  more  reliable  digital  filter  structures. 

During  the  contract  period  a number  of  technical  reports  on  the 
above  topics  were  written  and  sent  to  the  contract  monitor  at  WPAFB. 
These  reports  are  included  as  appendices,  ftany  of  the  results  described 
in  the  reports  are  briefly  mentioned  in  the  main  body  of  this  final 
report  and  the  reader  is  referred  to  the  appropriate  appendix  for 
further  information. 
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I.  DESIGN  OF  DIGITAL  FILTERS 


At  the  present  time,  the  majority  of  the  published  literature  on 
the  design  of  linear  time- invariant  digital  filters  describes  design 
using  frequency  domain.  The  most  probable  explanation  for  this  is  the 
wealth  of  knowledge  about  analog  filter  design.  From  the  frequency 
domain  point  of  view,  a desired  magnitude  function  is  typically  specified 
and  a discrete  transfer  function  is  found  whose  steady  state  magnitude 
function  satisfies  the  specifications.  However  it  is  also  possible  and 
in  many  cases  desirable  to  consider  time  domain  synthesis.  In  this 
section,  the  design  of  digital  filters  using  the  frequency  and  time 
domains  is  discussed. 

In  general  digital  filters  are  separated  into  two  classes, 
nonrecursive  or  finite  memory  and  recursive  or  infinite  memory  filters. 

A general  real  time  linear  time  Invariant  nonrecursive  digital  filter 
is  represented  by  the  difference  equation 
M 

Y[nT]  = y a,  X[(n  - k)T]  (1) 

k=l 

where  Y[nT]  is  the  filter  response  at  t = nT,  the  sequence 
is  the  sequence  of  filter  coefficients  and  the  X[ (n  - k)T]  are  input 
data  points.  In  (1)  T represents  the  sampling  Interval  which  is 
assumed  to  be  constant  and  M is  the  number  of  input  data  points  forming 
the  filter  window. 

The  general  expression  for  a real  time  linear  time  invariant 
recursive  digital  filter  of  order  N is  represented  by  the  difference 
equation 

M N 

Y[nT]  = I a,  X((n-k)T]-  [ b Y[(n-j)T].  (2) 

k=l  ^ j=l  J 
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Here  past  outputs  are  also  used  to  form  the  response  at  time  t = nT 
giving  rise  to  the  recursive  nature  of  the  filter.  The  design  task 
for  either  filter  is  to  determine  the  filter  coefficients  so  that  the 
filter  satisfies  specified  requirements. 

Design  Using  the  Time  Domain 

Synthesis  in  the  time  domain  is  based  on  a signal  model  whose 
output  is  typically  corrupted  by  noise.  For  nonrecurslve  filters  the 
most  common  signal  model  consists  of  a polynomial  over  the  finite 
window  which  implies  that  over  M data  points  the  signal  is  modeled 
by  a polynomial.  Because  this  type  of  digital  filter  is  useful  in  a 
variety  of  applications  such  as  radar  signal  processing,  the  general 
design  procedure  will  be  outlined.  First  the  concept  of  a transition 
matrix  for  a polynomial  signal  model  is  developed.  This  will  be  done 
only  for  a third  order  polynomial,  however,  the  results  are  easily 
generalized.  Let  a signal  X(t)  be  represented  by  a third  order 
polynomial  which  is  equivalent  to  the  differential  equation 


d X(t)  _ 


dt 


(3) 


If  a uniform  sampling  interval  is  assumed,  the  signal  and  the  first 
three  derivatives  at  the  point  t = (n  + h)T  can  be  written  as  Taylor 
series  expansions  about  the  point  t = nT  giving 

2 . . 3 . . . 

X[(n  + h)T]  -=  X[nT]  + hT  X[nT]  + X[nT]  + X [nT] , (4) 

X[(n  + h)T]  = X[nT]  + hT  X[nT]  + X [nT] , (5) 

X[(n  + h)T]  = X[nT]  + hT  X [nT]  (6) 

and 


j 


I 

i 

I 


4 


1 

, 

j 

■ 


2 
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X [ (n  + h)T]  = X [nT]. 

If  the  state  vector  is  defined  as 


(7) 


X[(n  + h)T] 


then  in  vector  form 


^X[(n  + h)T] 
X[(n  + h)T] 
X[(n  + h)T] 
\^X  [(n  + h)T]^ 


X[(n  + h)T]  = 1>[hT]  X[nT] 


where  ‘I>(nT]  is  the  state  transition  matrix 


$[hT]  = 


1 hT 


(hT)^  (hT) 
2!  3! 


3 ^ 


0 1 hT 


0 0 


0 0 


(hT)‘ 

2! 


hT 


0 


(8) 


(9) 


(10) 


To  eliminate  the  sampling  interval  T in  (10)  the  state  vector  is  now 
redefined  as 

^ X[(n  + h)!]"^ 

T X[ (n  + h)T] 

„2..  1 . (11) 

^ X[(n  + h)T] 

^ x'[(n  + h)T]  J 

Using  the  new  definition  the  vector  form  for  the  model  is 

X[(n  + h)T]  = $[h]  X[nT]  (12) 


X[(n  + h)T]  = 


where  the  new  state  transition  matrix  is  given  by 


i 


I 

I 


3 


2 3 

1 h h h 


0 1 2h  3h 


0 0 1 h 


0 0 0 1 


The  important  characteristic  of  $[h]  is  that 


4>[-h]  = [<I'[h]] 


The  determination  of  the  optimal  nonrecursive  filter  is  now 
formulated.  Given  that  the  signal  is  corrupted  by  noise,  the  observa- 
tion vector  at  time  t = nT  is  defined  by 


Y[nT]  = M X[nT]  + N[nT] 


where  Y[nT]  is  an  m element  observation  vector,  M is  a raxr  matrix 
that  designates  which  elements  of  the  state  vector  are  observed,  X[nT] 
is  an  r element  state  vector  and  N[nT]  is  an  m element  noise 
vector  with  zero  mean.  Given  Z observations,  the  total  observation 
vector  P[nT]  is  now  defined  as  a (m)(!l)  element  vector 


Y[nT] 


P[nT] 


Y[(n  - 1)T] 


Y[(n  - H - 1)T] 


The  design  problem  is  to  determine  the  weight  matrix  W that  makes  the 


estimate 


optimal . 


Z[nT]  = W P[nT]  . 
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In  order  to  write  ^[nT]  in  terms  of  the  state  vector  (15)  is 
substituted  into  (16)  giving 


P[nT]  = 


M X[nT] 


M X[ (n  - 1)T] 


M X[ (n  - i - 1)T] 


N[nT] 


N[(n  - 1)T] 


N[(n  - 5,  - 1)T] 


Substitution  of  (12)  into  (18)  now  allows  P[nTl  to  be  written  as 


P[nT]  = T X[nT]  + u[nT] 


where  u[nT]  is  the  total  noise  vector  on  the  right  side  of  (18).  In 
(18)  the  matrix  T is  given  by 


M 4>[h] 


M 4>[h  -8,-1] 


and  as  will  be  shown  later  the  T matrix  plays  an  essential  role  in 
the  filter  design. 

Finally  the  optimal  estimate  is  assumed  to  be  an  unbiased  estimate 
which  results  in  a constraint  relationship  between  W and  T.  To 
derive  this  (19)  is  substituted  into  (17)  giving 

Z[nT]  = WT  X[nT]  + W u[nT]  (21) 

Since  ^[nT]  is  an  unbiased  estimate  of  X[nT],  taking  the  mean  of  (21) 
now  yields  the  equation 


X[nT]  = WT  X[nT] 


which  is  only  satisfied  if  the  constraint 


WT  = I 


is  valid. 


5 


' For  a conventional  least  squares  estimate  the  minimization  problem  is 

handled  by  minimizing  the  cost  function 

E[nT]  = [T  Z[nT]  - P[nT]  ] [T  Z[nT]  -P[nT]]  (2A) 

Note  that  if  the  total  noise  vector  is  zero,  the  cost  function  is  zero 
i and  any  weight  matrix  W that  satisfies  (23)  will  result  in  an  optimal . filter . 

i 

■ Since  ^[nX]  is  the  unknown,  the  differentiation  of  (24)  with  respect 

to  ^[nl]  and  the  setting  of  the  result  to  zero  gives 

I'^T  Z[nT]  - P[nT]  = 0.  (25) 

A 

The  solution  of  (25)  gives  an  optimal  estimate  denoted  by  X[nX] 
which  is 

X[nl]  = [t‘^T]“^  P[nT]  (26) 


so  that  the  optimal  weight  vector  is 


“ ^ — 1 f 

W = [XT]  X 


In  practice,  the  X matrix  can  often  be  written  by  inspection. 


For  example  since  usually  only  the  data  and  not  the  derivatives  are 


known,  the  observation  is  a scalar.  For  this  important  case,  the  X 
matrix  for  a three  point  window  and  a second  order  polynomial  signal  is 


given  by 


Since 

only  tlie 

data  X 

[nT] 

is 

available  M is  defined  as  the  ma 

tr  ix 

11  0], 

. The  r 

matrix 

is  now  given  as 

j 

1 0 

- 

i 

[1  0] 

1 

-1 

T = 

0 

1 

(30) 

[1  0] 

1 

-2 

; 

■ 

0 

1 _ 

i 

! which 

when  multiplied  out  gives 

(28).  The  optimal  weight  matrix 

is  now 

given 

by 

1 

5 2 

-1 

"-t 

3 0 

-3, 

(31) 

Using  (31)  the  two  elements  of  the  optimal  estimate  vector  are  estimates 
of  the  data  and  the  first  derivative  and  are  given  by 

X^[nT]  = I Y[nT]  + | Y[ (n  - 1)T]  - | Y[ (n  - 2)T]  (32) 

and 

X2[nT]  = Y[nT]  + 0 Y[ (n  - 1)T]  - ^ Y[(n  - 2)T].  (33) 

For  this  formulation,  the  T matrix  can  be  easily  generalized. 

The  number  of  rows  in  the  T matrix  is  the  size  of  the  data  window. 

The  first  column  consists  of  ones,  and  the  elements  in  all  other  columns 
are  given  by 

T_  = (-1)^  (i)^  j ,i  1 0 (34) 

Thus  for  a five  point  filter  the  T matrix  for  a third  order  polynomial 
fit  has  five  rows  and  three  columns  giving 

f 

I 
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Given  the  optimal  estimate  formed  using  conventional  least  squares, 
the  covariance  matrix  for  the  estimate  is 

C = W*^RW  (36) 

where  R is  the  covariance  matrix  of  the  total  noise  vector  u[nT]. 

In  most  applications  the  noise  is  assumed  to  be  statistically  time 
invariant  so  the  elements  of  R are  all  constant.  Carefully  note  that 
conventional  least  squares  does  not  consider  any  correlation  between 
noise  samples  and  the  estimate  is  not  optimum  if  correlation  exists. 
However  in  many  cases  the  estimate  is  satisfactory  and  the  optimal 
weight  matrix  is  easily  found  via  a computer.  In  tenns  of  a good  design 
tool,  (27)  is  easily  programmed  and  with  an  Interactive  graphics  terminal 
the  designer  can  quickly  determine  filter  coefficients  by  simply  giving 
the  number  of  row  and  column  elements  of  the  T matrix. 

If  the  measurement  noise  is  correlated  from  sample  to  sample,  the 
optimal  estimate  is  found  using  the  concept  of  minimum  variance  or 
weighted  least  squares.  The  resulting  optimal  weight  matrix  for  the 
estimate  is  given  by 

W[mln  var]  = R~^T]”^  R~^  (37) 

where  R is  the  autocovariance  matrix  of  the  measurement  noise.  The 
corresponding  covariance  matrix  for  the  estimate  is  given  by  (36)  using 
the  weight  matrix  of  (37). 
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The  extension  of  the  theory  to  time  invariant  recursive  filters  is 
not  simple  because  the  covariance  matrix  of  the  estimate  vector  is  time 
varying  until  steady  state  is  reached.  As  a result  little  information 
in  the  literature  is  available,  however,  the  problem  has  been  investigated 
by  Brubaker  and  Harper  [1]  for  first  and  second  order  recursive  filters. 

For  the  case  where  the  state  equation  is  driven  by  input  noise  the  well 
known  Kalman  filter  results.  This  filter  has  been  the  subject  of 
numerous  papers,  and  no  theory  will  be  given  here.  However  its  implemen- 
tation is  not  as  well  known  and  the  report  Lhau  describes  a modified 
Kalman  filter  where  the  effects  of  product  rounding  errors  and  random  errors 
in  the  filter  coefficients  are  considered  is  being  sent  under  a separate  cover. 

Because  time  domain  synthesis  methods  have  great  potential  in  a 
variety  of  applications  of  interest  to  the  Air  Force,  it  is  anticipated 
that  new  procedures  for  designing  time  invariant  recursive  filters  will 
be  developed  in  the  future.  The  advantage  to  be  gained  is  greater 
noise  reduction  with  less  computing.  A good  reference  on  time  domain 
synthesis  is  by  Morrison  [2]  who  gives  a good  background  in  current 
time  domain  s3mthesis  procedures. 

Design  Using  the  Frequency  Domain 

The  design  of  nonrecursive  filters  via  the  frequency  domain  normally 
starts  with  a specified  magnitude  function  for  the  digital  filter.  The 
task  is  to  determine  the  coefficients  in  (1)  so  that  the  magnitude 
function  for  the  filter  satisfies  the  specifications.  The  magnitude 
function  is  found  by  first  taking  the  Z transform  of  (1)  giving 

M -k 

Y(Z)  = I a,^  Z X(Z)  . (38) 

k=l 

The  discrete  transfer  function  is  now  given  by 
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J 


H(Z)  = I \ (39) 

k=r 

and  the  steady  state  magnitude  function  is  found  by  letting  Z = 
and  taking  the  absolute  value  of  both  sides  of  (39).  This  yields 

k=l 

One  method  of  determining  the  coefficients  utilizes  linear 
programming  and  is  described  in  references  [3,  4].  Another  procedure 
described  by  Farden  and  Scharf  [5J  uses  a Wiener  filter  structure  to 
minimize  the  noise  while  giving  a good  approximation  to  the  magnitude 
function.  It  is  also  possible  to  use  a truncated  Fourier  expansion  of 
the  periodic  magnitude  function  since  the  frequency  function  of  any 
discrete  filter  is  periodic  in  the  frequency  domain. 

The  primary  problem  with  the  above  methods  for  design  is  the  large 
size  of  the  data  window  that  is  required  to  synthesize  any  filter  with 
reasonably  sharp  attenuation  characteristics.  For  example  to  obtain 
sharp  attenuation  between  the  pass  and  stop  bands  in  reference  [5] 
window  sizes  of  128  to  256  coefficients  were  required.  The  resulting 
filters  require  too  much  computational  time  to  be  useful  in  medium  to 
high  frequency  applications.  As  a result  the  nonrecurslve  filter 
appears  to  most  suitable  for  time  domain  synthesis  where  fewer  coeffi- 
cients are  utilized  and  the  filter  is  used  for  preprocessing. 

For  recursive  digital  filters  there  are  a variety  of  frequency- 
domain  procedures.  First  Rabiner  [6]  has  extended  the  linear  proj  ramming 
procedure  for  nonrecursive  filters.  However  the  procedure  does  not 
appear  to  be  satisfactory  if  the  filter  has  sharp  attenuation  between 
the  pass  and  stop  bands.  Also  trigonometric  polynomials  are  employed 
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in  a procedure  similar  to  the  design  of  Butterworth  and  Chebyshev 
filters  in  the  analog  domain  [7],  Other  methods  use  impulse  invariance 
and  frequency  sampling  [7]. 

In  general  none  of  the  above  methods  has  gained  wide  acceptance. 
One  reason  for  this  is  the  lack  of  phase  information  which  is  important 
in  digital  control  work.  However  the  primary  reason  deals  with  the 
implementation.  When  high-order  recursive  filters  are  implemented 
directly,  the  accuracy  with  which  the  coefficients  must  be  represented 
increases  rapidly  as  the  filter  order  increases.  The  same  effect  is 
present  when  implementing  high-order  analog  filters  and  a general 
design  principle  is  that  analog  filters  should  be  implemented  as  a 
series  or  parallel  connection  of  first  and  second  order  sections.  The 
same  implementation  strategy  holds  for  digital  filters,  however,  with 
high  order  filters  the  polynomial  factoring  must  be  done  with  great 
accuracy. 

The  use  of  the  bilinear  Z transform  on  the  other  hand,  allows 
the  implementation  directly  from  the  series  or  parallel  connection  of 
the  corresponding  analog  filter  since  the  bilinear  Z transform  of  a 
sum  or  product  is  the  sum  or  product  of  bilinear  Z transforms.  Note 
that  this  is  not  true  of  the  classical  Z transform.  For  example, 
if  the  Z transform  of  two  second  order  analog  filters  in  cascade  is 
to  be  determined,  the  sections  must  first  be  multiplied  to  give  a 
fourtli  order  transfer  function.  Then  the  weighting  function  is  found 
by  taking  tlie  inverse  Laplace  transform.  Finally  the  Z transform  of 
the  weighting  function  is  computed.  In  general  this  is  not  a pleasant 
task. 

Because  the  bilinear  Z transform  method  is  widely  used  and  at 
present  is  probably  most  appropriate  for  most  Air  Force  applications 
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this  is  the  only  frequency  design  method  for  recursive  filters  that  is 
considered  in  this  report.  However  all  of  the  results  concerning  the 
implementation  apply  directly  to  filters  designed  using  other  methods. 
Note  that  in  digital  control,  design  using  the  bilinear  Z transform  is 
called  lustin' s method. 

The  bilinear  Z transform  method  first  requires  the  design  of  an 
analog  filter  that  meets  the  specifications.  If  the  transfer  function 


for  the  analog  filter  is  represented  as  H(s)  the  discrete  transfer 
function  for  the  corresponding  digital  filter  is  found  by  the  substitution 


2 Z - 1 ,,, 
^ ‘ T TVl 

where  T represents  the  sampling  interval.  In  many  references  (41) 
is  called  the  extended  bilinear  Z transform.  The  procedure  is  now 
illustrated  by  a very  simple  example. 


Consider  the  first  order  analog  filter  with  a dc  gain  of  one  and 
a transfer  function 


H(s)  = 


s + a 


Using  the  bilinear  Z transform  or  Tustin's  method,  (41)  is  substituted 


into  (42)  to  give  the  discrete  transfer  function 


H(Z)  = 


a (Z  + 1) 
o 

Z - b. 


where  the  constants  a^  and  b^^  are  given  by 


o 2 . 

+ a 
T 


, T " ® 
1 ■ ^ 
Y + a 
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In  terms  of  the  implementation,  as  the  sampling  time  becomes  small,  it 
obviously  takes  more  decimal  digits  or  binary  bits  to  accurately  represent 
Che  coefficients.  This  effect  is  magnified  in  higher  order  digital 
filters  so  that  a very  important  design  consideration  is  to  choose  the 
lowest  possible  sampling  rate  which  is  equivalent  to  the  largest  possible 
sampling  interval. 

It  appears  that  the  use  of  the  bilinear  Z transform  that  does  not 
contain  the  2/T  factor  might  eliminate  the  dependence  of  the  coeffi- 
cients on  the  sampling  interval  T,  however  the  warping  effect  cancels 
out  any  benefits.  This  warping  using  the  bilinear  Z transform  of  (41) 
is  now  illustrated.  For  a given  digital  frequency  the  corresponding 

analog  frequency  is  represented  by  u . In  steady  state  Z = 

3i 

and  s = jw  so  substitution  into  (41)  and  manipulating  gives 

Wg  = -|  Tan(o)^  -|)  (46) 

Thus  any  frequency  lo^  for  the  digital  filter  has  a corresponding 
analog  frequency  that  is  warped  by  the  relationship  of  (46).  This 

implies,  for  example,  that  the  critical  frequencies  such  as  the  -3db 

frequencies  are  not  the  same  in  the  analog  and  digital  filters.  Note 
that  if  the  2/T  factor  is  not  used,  the  warping  will  be  much  more 
severe. 

To  minimize  this  warping  the  sampling  interval  T can  be  chosen 
small  enough  so  that  the  approximation 

Tan  “d  ^ 

holds.  This,  however,  increases  the  sampling  rate  for  the  digital 
filter  which  in  turn  requires  greater  accuracy  for  coefficient 
representation.  In  terms  of  the  computer  a longer  word  length  is 
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needed  for  coefficient  representation.  Another  method  to  avoid  warping 


at  critical  frequencies  is  to  prewarp  the  critical  frequencies  of  the 
analog  filter  before  the  bilinear  Z transform  is  applied.  However, 
other  frequencies  in  the  filter  are  still  warped  with  respect  to  the 
analog  and  digital  frequencies. 

Two  design  programs  utilizing  the  bilinear  Z transform  have  been 
sent  to  WPAFB  as  part  of  the  current  work.  The  first  package  is  used 
for  the  design  of  digital  Butterworth  and  Chebyshev  band-pass  and  band- 
stop  filters  with  up  to  six  second  order  sections  in  cascade.  The  user 
need  only  enter  the  type  of  filter,  the  order,  the  sampling  rate  and  the 
upper  and  lower  -3db  frequencies.  A description  of  the  program  is 
given  in  Appendix  A. 

The  second  program  is  used  for  designing  Butterworth  and  Chebyshev 
low  pass  digital  filters.  Here  the  designer  enters  the  type  of  filter,  j 

the  filter  order,  the  sampling  rate  and  the  -3db  frequency.  The 
output  consists  of  the  coefficients  for  each  second  order  section  in  the 
cascade.  A description  of  the  program  is  given  in  Appendix  B. 

Sampling  Rate  Selection 

One  of  the  key  design  parameters  in  digital  filter  design  that  is 
usually  overlooked  is  the  sampling  rate.  In  many  designs  the  sampling 

rate  is  chosen  higher  than  necessary  which  can  lead  to  a variety  of  ^ 

problems  such  as  instability  due  to  coefficient  rounding.  One  argument 
for  a high  sampling  rate  is  the  need  for  analog  filtering  before  sampling. 

This  is  used  to  prevent  aliasing  of  high  frequency  noise.  If  the  analog  \ 

filter  is  to  cause  minimal  phase  shift  in  the  signal  band  and  good  | 

attenuation  at  the  sampling  frequency  it  is  necessary  to  make  the  I 

sampling  frequency  high  with  respect  to  the  maximum  signal  frequency.  1 


lA 


Although  the  magnitude  characteristic  for  any  digital  filter  is 


periodic  in  frequency,  random  noise  can  still  be  reduced  although  the 

amount  of  reduction  is  less  than  that  of  an  equivalent  analog  filter. 

To  provide  some  insight  into  the  noise  problem  and  its  relationship 

to  sampling,  consider  the  filter  given  by  (42)  with  a = 1.  For  a noise 

2 

input  with  variance  a and  an  autocorrelation  function 


R (t)  = a e 

XX 


2 


(48) 


the  input  power  spectral  density  is 
0^  2(0 

"xx(“)  = -T--T 

(Oj^  + (O 


(49) 


where  to^^  is  the  -3db  bandwidth  of  the  noise  whose  dc  value  is 
2 

2o  /(Oj^.  The  variance  of  the  output  noise  for  the  analog  filter  is 

2 


Var  {Y}  = - - - 

a (0,  + 1 


(50) 


Note  that  if  co^  = 1 the  filter  has  reduced  the  variance  of  the  noise 
by  one-half.  If  the  noise  is  sampled  to  drive  a corresponding  digital 
filter  designed  using  the  bilinear  Z transform,  the  variance  of  the 
noise  at  the  digital  filter  output  is 


Var.{Y^}  = - 
d n 2 


ahl 


■“iT) 


Y(l-e 


-(0]^T 


) + (1  + e 


-(OiT. 


(51) 


As  T approaches  zero  the  two  variances  given  by  (50)  and  (51)  converge. 

However  for  a given  sampling  rate  the  digital  filter  will  reduce  the 

variance  of  the  noise  for  any  noise  bandwidth  co^^.  This  is  shown  in 

Table  1 for  a noise  bandwidth  of  ten  radians  per  second  which  Is  ten 

times  the  filter  bandwidth.  In  the  table  is  the  sampling  frequency 

defined  as  (o  = 2tt/T. 
s 
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Table  1 

EVALUATION  OF  EQUATIONS  50  AND  51  FOR  A NOISE  BANDWIDTH  OF 
TEN  RADIANS  PER  SECOND 


T 

CJ 

s 

Var  (Y) 
a 

a n 

2 

2 

0 

0 

2.0 

IT 

0.091 

0.500 

1.0 

2tt 

0.091 

0.333 

0.5 

4tt 

0.091 

0.202 

0.25 

8tt 

0.091 

0.128 

Looking  at  the 

table , 

the 

difference  between 

the 

noise  levels  of 

the  analog  and  digital  filters 

at  a sampling  rate 

of 

4tt  times  the 

filter  bandwidth  is 

about 

6db, 

, In  general,  unless  a 

serious  noise 

problem  is  present,  no  analog  prefiltering  may  be  needed.  However  if 
filtering  is  done,  wide-band  filtering  should  be  tried  or  a linear  phase 
filter  should  be  used  to  allow  digital  phase  compensation  to  be  employed. 
This  can  be  done  via  the  use  of  predictive  nonrecursive  filtering. 


II.  IMPLEMENTATION  OF  DIGITAL  FILTERS 


In  this  section,  the  procecures  that  can  be  used  by  the  designer 
for  determining  the  required  computer  word  length  needed  for  the 
implementation  of  a digital  filter  are  discussed.  This  topic  has  been 
well  developed  in  reports  included  as  appendices  C,  D,  and  E and  the 
reader  should  refer  to  these  appendices  for  more  infomation. 

In  a digital  filter  error  is  introduced  by  rounding.  For  fixed 
point  aritlimetic,  errors  are  caused  by  the  rounding  of  input  data,  the 
filter  coefficients,  and  the  intermediate  products.  For  floating  point 
arithmetic,  additional  error  is  caused  by  the  round*ng  of  sums.  Because 
fixed  point  arithmetic  is  inexpensive,  fast  and  reliable  it  is  the  most 
suitable  for  aircraft  systems.  However,  when  fixed  point  arithmetic  is 
used,  the  designer  must  carefully  specify  the  amount  of  rounding  error 
that  is  acceptable  and  then  must  determine  the  configuration  and  the 
word  length  for  the  filter  that  will  allow  the  specifications  to  be 
met.  In  practice  this  must  also  be  done  for  floating  point,  however, 
most  designers  do  not  undertake  the  task.  As  a result,  the  specified 
computer  requirements  are  almost  always  excessive  resulting  in  more 
expensive  less  reliable  aircraft  computer  systems. 

Because  fixed  point  arithmetic  is  satisfactory  for  most  aircraft 
applications,  this  arithmetic  will  be  considered  in  this  report. 

However,  the  results  can  easily  be  extended  to  floating  point  to  find 
the  exact  word  length  requirements  for  a floating  point  system. 

Rounding  of  the  Input  Data 

Any  analog-to-digital  converter  has  an  output  that  is  represented 
in  the  computer  by  a finite  number  of  bits.  Since  most  analog  inputs  to 
a converter  are  bipolar,  most  analog-to-digital  converters  are  specified 


to  have  a two's  complement  binary  output  consisting  of  P bits  and  a 
sign  bit.  This  means  that  each  data  point  is  rounded  with  an  error 


The  binary  equivalent  of  q is  1.  If  the  output  of  the  A/D  converter 

is  scaled  to  a binary  fraction,  the  value  of  q in  signal  units  is 

still  10  millivolts,  however,  the  binary  equivalent  is  given  by  2 

Input  rounding  or  quantization  has  been  investigated  by  a number 

of  researchers.  Bennett  [8]  and  Widrow  [9]  have  shown  that  for  a rounding 

form  of  quantization,  the  errors  can  be  modeled  as  independent  random 

2 

variables  with  zero  mean  and  variance  q /12.  This  model  has  been 
experimentally  verified  for  input  signals  that  traverse  several 
quantization  levels  from  sample  to  sample.  The  truncation  form  of 
quantization  is  similar  except  that  the  mean  value  for  each  error  is 
q/2.  Note  that  in  terms  of  probability  the  rounding  error  is  assumed 
to  be  unifonnly  distributed  about  the  mean. 

For  a random  model,  input  rounding  consists  of  an  additive  noise 
source  at  the  input.  The  variance  of  the  error  at  the  filter  output 
is  given  by 

2 2 

Var{Y  } = ^ y H [kT]  (54) 

k=0 

where  H[kT]  is  the  inverse  Z transform  of  the  transfer  function  for 
the  digital  filter.  Note  that  (54)  is  time  varying  and  the  variance 
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becomes  a constant  in  steady  state  by  letting  n approach  infinity. 

In  this  case  the  variance  is  also  given  by  the  discrete  Wiener-Khinchine 
relationship 

2 

Var{Y[==]}  = ^ ^ ^ H[Z]  z"  dZ  (55) 

Given  a specified  value  for  the  variance  at  the  filter  output  due  to 
input  rounding,  the  word-length  of  the  filter  can  be  determined  from 
(54)  and  (55).  This  is  done  as  follows. 

Since  the  discrete  transfer  function  is  the  ratio  of  two  polynomials 
in  Z,  the  inverse  Z transform  of  the  transfer  function,  called  the 
weighting  sequence  is  found  via  polynomial  division.  By  squaring  each 
term  and  adding, the  asymptotic  steady  state  variance  is  found.  The 
steady  state  variance  can  also  be  found  using  the  procedure  developed  in 
Appendix  3.  This  formulation  is  easily  programmed  for  use  with  an 
interactive  graphics  terminal. 

Given  the  summation  and  a specified  value  for  the  variance,  q is 
given  by 

12  (Specified  Variance) 

q < 5 — t 

^ — CO 

y H^[kT] 

L k^O 

Since  the  specified  variance  is  usually  given  as  a signal  to  noise 
ratio,  tlie  signal  range  is  part  of  the  specified  variance  and  the  actual 
number  of  bits  can  be  determined.  For  example  if  a sinusoidal  waveform 
with  maximum  amplitude  is  defined  as  the  signal  output,  then  the  maximum 
signal  power  in  an  analog  sense  is 

SP  = — . (57) 

/2 


1/2 


(56) 
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The  noise  power  is 


I 2 


and  for  a s/N  ratio  of  say  -20  db 


-20  _<  10  log 


q(2^-l) 

/2 

I 1 - ^ 

h 5;  H [kT] 

k-o 


where  the  range  R is  (2^-l)q  and  p is  the  number  of  binary  bits. 
Solving  (56)  for  ^ now  yields 


2^  J (4.08)10"^/  S H^fkT]  + 1 


For  a given  weighting  sequence,  p is  found  to  be  the  smallest  number 
that  satisfies  the  inequality. 

It  is  also  possible  to  determine  q and/or  p via  use  of  a 
Gaussian  distribution  for  the  noise.  Here,  the  designer  assigns  a 
confidence  level  for  an  acceptable  error.  Use  of  standardized  tables 
then  allows  the  value  of  p to  be  determined.  This  procedure  is 
outlined  in  more  detail  in  Appendix  D. 

If  a statistical  analysis  is  not  satisfactory,  then  worst  case 
design  procedures  can  be  Initiated.  A method  for  finding  the  maximum 
upper  bound  on  the  error  due  to  input  rounding  was  Introduced  by 
Bertram  [10].  In  this  analysis,  the  error  is  represented  as  the 
convolutional  sum 


e[nT]  = Z H[kT]  N [(n-k)T] 
k=0 
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where  the  input  noise  samples  are  defined  as  the  sequence  {N[kT]!. 
Taking  the  absolute  value  of  both  sides  and  choosing  the  maximum 
absolute  value  of  the  error  gives  the  maximum  upper  bound  for  the  error 
at  the  filter  output  as 


|e[nTl  \ r.  lH[kTl  | = 


An  alternative  procedure  has  been  developed  by  Slaughter  [11]. 
Here  each  error  source  is  considered  as  a step  function  with  height 


q/2.  The  error  is  then  given  as  the  dc  steady  state  output.  In 
practice,  the  determination  of  the  word  length  using  a confidence  level 
of  95  percent  usually  results  in  a smaller  word  length  requirement.  The 
method  of  Slaughter  gives  a smaller  word  length  than  that  of  Bertram's. 

In  most  applications,  simulation  shows  the  Bertran's  method  yields 
pessimistic  results.  Examples  for  each  procedure  are  given  in  Appendix  D. 
with  the  results  including  error  due  to  product  rounding. 


Rounding  of  Products 

In  any  digital  filter,  the  multiplication  of  a P bit  and  sign 
number  by  an  m bit  and  sign  coefficient  results  in  a product  represented 
by  (P  + m)  bits  and  sign.  In  many  systems,  p and  m are  assumed  to 
be  equal  so  the  product  contains  2p  bits  plus  sign.  In  most 
implementations,  the  products  are  rounded  resulting  in  an  error  e 


bounded  by 


e < ^ 

P - 2 


where  q is  the  least  significant  bit  of  the  rounded  product.  In 


practice  this  value  of  q^  may  not  be  the  same  for  various  products  and 
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not  the  same  as  the  q for  the  input  rounding  error.  However,  for 


the  purpose  of  analysis  all  values  of  q will  be  taken  as  the  same. 

As  is  the  case  with  input  rounding  error,  product  rounding  errors 

are  assumed  to  be  uncorrelated  zero-mean  random  variables  with  variances 
2 

q /12.  Each  product  rounding  error  is  assumed  to  be  an  additive  noise 
source  following  the  multiplication  in  the  filter  block  diagram.  The 
only  difference  in  the  analysis  is  the  fact  that  the  product  rounding 
errors  occur  Inside  of  the  filter.  As  a result,  the  form  of  implementa- 
tion plays  an  important  factor  in  the  implementation  design  phase. 

Currently,  most  digital  filters  are  implemented  using 
the  Direct  Form  1 (DFl)  or  Direct  Form  2 (DF2) . Block  diagrams  for  a 
third  order  filter  are  shown  in  Figs.  3 and  4 of  Appendix  D.  In 
addition,  it  is  well  known  that  to  minimize  the  effect  of  product 
rounding,  a high  order  filter  should  be  Implemented  as  a cascade  or 
parallel  connection  of  first  and  second  order  sections.  This  is 
pointed  out  by  Liu  [12]  and  in  Appendix  C.  For  example  in  Appendix  C. 
the  autocovariance  sequence  for  a fourth  order  Butterworth  filter 

2 

implemented  in  DFl  is  shown  in  Fig.  2.  The  variance  is  about  48.75q  . 

Note  that  the  variance  is  given  when  k-=0.  For  two  second  order 

2 

sections  in  cascade  the  variance  is  5.83q  whlcli  is  a vast  improvement. 

Another  problem  that  is  present  and  closely  related  to  product 
rounding  is  overflow.  While  a certain  implementation  structure  may 
result  in  minimal  product  rounding,  when  the  necessary  scaling  is 
employed  to  prevent  overflow  the  structure  can  exhibit  increased  error 
due  to  product  rounding.  In  practice  the  overflow  often  occurs 
during  the  transient  response  of  the  filter  to  a given  input  and  it  is 
difficult  by  analytical  means  to  determine  when  overflow  is  present. 
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As  a result,  the  design  of  a good  filter  structure  that  is  satisfactory 
from  both  product  rounding  and  overflow  almost  always  requires  the  use 
of  interactive  design  software  that  allows  the  designer  to  simulate 
various  configurations  over  the  specified  set  of  input  signals. 

To  demonstrate  the  ideas  in  a single  example,  a first  order 
digital  filter  will  be  implemented  in  DFl  and  DF2  as  shown  in  Figs,  la 
and  lb.  The  transfer  function  for  the  filter  is  given  by 


H(Z) 


"^0^'  ^1  ^ N(Z) 

Z - b^  D(Z) 


(63) 


Referring  to  Fig.  la,  the  three  multiplication  errors  are  operated  on 
only  by  the  pole  so  tlie  error  is  given  by  the  convolutional  sum 


eDFitnT]  = I b^  ic^fCn-k)!]  + c^((n-k)T]  + 6^((n-k)T]}  (64) 
k=0 


In  (62)  b^  is  the  Inverse  Z transform  at  1/P(Z)  in  (61).  The 


steady  state  variance  is  now  given  by 


Var  {epp^[”l} 


- 


(65) 


In  terms  of  overflow,  for  a unit  step  input  the  response  in  the 
domain  is 


y(2) 


^0^'  ^1 
(Z-b^) 


z 

Z-1  ■ 


z 


(66) 


The  corresponding  time  response  in  steady  state  is 

a + a 

(y  ) =-^ 


(67) 


Thus  for  unity  gain  (65)  is  set  equal  to  one  and  for  this  case  no 
overflow  will  occur  in  the  DFl  implementation. 


23 


For  the  Direct  Form  2,  the  variance  of  the  error  In  steady  state  is 


„ 2 1 (a  2+2a  a b +a 

Var{eDF2fnT]l  = -^  + 7—2 


(68) 


1-b 


1 


For  a dc  gain  of  one  substitution  of  this  into  (66) 

1-b^  ° ^ ^ 

for  a^  = a^  = — gives 


^ 2 2 ri-b  ) 

Var{eDF2[nT]}  = ^ + “I 


(69) 


Obviously  without  scaling  (67)  is  less  than  (63)  for  values  of  b 


1 


between  -1  and  1.  However,  for  a step  input,  the  value  of  T for 

n 

the  DF2  in  steady  state  is 


(T  ) 


n n->«'  1-b, 


(70) 


For  positive  values  of  b^^,  (68)  is  greater  than  one  and  overflow 
will  occur.  To  avoid  this,  the  input  is  usually  scaled  down  by 

a scale  factor  l-b^^.  This  additional  multiplication  introduces 
another  noise  term  at  the  input.  To  obtain  the  correct  scaling  at  the 
output  the  two  coefficients  a^  = aQ/(l-b^)  and  aj^  = aj^(l-bj^). 

For  ^0  ~ ^1  ~ (l-bj^)/2,  the  resulting  filter  has  a steady  state 
variance 


2 2 2 2 

Var{e^„.(nT) } = ^ (T^^“)  + (tT^) 


'DF2 


12  ■ 12  ' 12  '1-b^ 


(71) 


where  r is  the  scale  factor  for  the  quantization  caused  by  the  scale 
factor  at  the  input.  In  general  r is  greater  than  one.  A comparison 
of  the  variance  of  the  two  forms  in  the  scaling  as  employed  in  the  DF2 
implementation  is  given  in  Table  2.  A dc  gain  of  one  is  assumed  and 

■■■0  ■ ",  ■ 
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TABLE  2 


COMPARISON  OF  THE  VARIANCE  OF  DFI  AND  DFZ 


”l 

^0 

r 

Var  DFI 

7 

q /12 

Var  DF2 
q“/12 

. 1 

.45 

.45 

2 

3.03 

4.77 

. 2 

.4 

.4 

2 

3.13 

5.13 

. 3 

.35 

.35 

2 

3.30 

5.57 

.4 

.3 

.3 

2 

3.57 

6.17 

.5 

.25 

.25 

4 

4.0 

19.00 

.6 

_ 2 
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4.69 

23.25 

.7 

.15 

. 15 

4 

5.88 

30.33 

.8 

. 1 

. 1 

8 

8.33 

164.00 

.9 

.05 

.05 

16 

15.79 

1287.00 

.99 

.005 

.005 

128 

150.75 

819252.00 

From  the  table  it  is  obvious  that  when  scaling  is  employed  in  the 
first  order  filter,  the  output  error  due  to  product  rounding  is  much 
greater  for  the  DF2  form.  In  general,  scaling  must  be  done  via  the  use 
of  interactive  graphics  for  the  desired  set  of  input  signals.  Also 
it  is  important  to  recognize  that  for  a given  filter  there  is  no 
analytical  method  to  determine  the  best  form  for  a given  section  and 
an  optimal  ordering  of  the  sections  that  will  result  in  no  overflow 
and  minimum  variance  due  to  product  rounding.  Thus,  the  designer  must 
carefully  utilize  interactive  computation  along  with  a graptiics  terminal 
to  scale  the  filter  and  order  the  sections. 

More  detail  on  the  implementation  of  digital  filters  is  given  in 
Appendix  . While  the  procedures  described  in  the  Appendix  are 
pointed  toward  digital  control,  they  are  directly  applicable  to  any 
digital  filter.  Both  input  and  product  rounding  errors  along  with 
scaling  are  discussed.  The  method  for  computing  the  variance  is  shown 
in  Appendix  along  with  an  example  for  product  rounding.  This  algorithm 
is  well  suited  for  use  with  an  interactive  graphics  terminal. 

One  other  phenomenon  that  is  present  due  to  product  rounding  is 
the  limit  cycle.  At  present  there  is  no  general  theory;  however. 
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it  is  well  known  that  rounding  or  truncation  of  products  can  cause  a 
digital  filter  to  exhibit  an  oscillation  or  limit  cycle.  This  occurs 
with  both  fixed  and  floating  point  arithmetic.  A reference  is  by 
Brubaker  and  Gowdy  |13]. 

CoefficienL  Rounding 

The  last  source  of  error  in  a filter  is  caused  by  the  rounding  or 
truncation  of  the  filter  coefficients.  While  some  writers  describe 
these  errors  as  random,  in  practice  they  are  deterministic  and  under 
control  of  the  programmer.  In  filters  with  rigid  specifications,  the 
actual  rounding  or  trucation  of  the  coefficients  must  be  done  with 
care  or  the  filter  charcteristics  will  not  satisfy  the  specifications. 

For  example,  for  low-pass  filters  the  dc  gain  can  change  dramatically 
if  the  coefficients  are  not  handled  properly.  This  is  illustrated  in 
Appendix  E where  a second  order  lag-lead  digital  filter  is  Investigated. 
By  use  of  a strategy  for  rounding  and  truncating  the  coefficients,  the 
filter  coefficients  can  be  handled  with  twelve  bits  plus  a sign  bit. 

When  conventional  rounding  is  employed  the  use  of  twelve  bits  plus 
sign  gives  a dc  gain  of  4.88  x 10^  when  the  desired  dc  gain  is  one. 

At  present  there  are  not  many  general  procedures  for  coefficient 
rounding.  In  Appendix  E,  the  rounding  is  handled  by  first  specifying 
how  much  error  in  the  filter  magnitude  and  phase  functions  is  acceptable. 
Then,  the  use  of  a differential  approximation  allows  the  error  to  be 
represented  in  terms  of  the  filter  coefficients  and  the  coefficient 
errors.  A linear  programming  method  is  then  used  to  establish  a region 
where  the  coefficient  error  vector  must  lie.  Note  that  while  the 
strategy  in  Appendix  E is  used  for  digital  control  algorithms  it  is 
applicable  to  any  filter. 
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Another  important  consideration  that  was  pointed  out  in  the  design 
section  is  the  fact  that  the  accuracy  with  which  the  coefficients  must 
be  represented  is  dependent  on  the  sampling  interval  T.  For  high 
order  filters  the  accuracy  requirement  drastically  Increases  as  the 
order  of  the  filter  increases.  As  a result,  another  reason  for 
Implementing  filters  as  a cascade  or  parallel  connection  of  first  and 
second  order  sections  is  to  reduce  error  due  to  coefficient  rounding. 
Thus  the  designer  must  carefully  compromise  the  sampling  time  to  meet 
the  specifications  and  to  minimize  the  effect  of  coefficient  rounding. 
Then  the  filter  is  implemented  as  a cascade  or  parallel  connection  of 
first  and  second  order  sections.  The  actual  determination  of  the 
number  of  bits  needed  for  each  first  and  second  order  section  is  done 
via  the  strategy  described  in  Appendix  E. 

To  summarize  the  overall  implementation  strategy,  the  designer 
first  designs  a given  filter  using  the  lowest  possible  sampling 
frequency.  The  filter  is  then  factored  for  implementation  as  a 
cascade  or  parallel  connection  of  first  and  second  order  sections. 

The  ordering  of  the  sections  and  the  type  of  connection  to  be  utilized 
is  done  using  interactive  computer  packages  to  establish  a configuration 
that  minimizes  the  effect  of  product  rounding  and  still  has  no  overflow. 
Then  the  error  due  to  input  rounding  and  product  rounding  is  used  to 
compute  the  word  length  needed  to  meet  the  specifications  on  these  two 
errors.  Finally  the  coefficients  are  rounded  and  tlie  word  lengtti  is 
determined  to  meet  the  specifications  on  coefficient  rounding.  The 
final  word  length  is  chosen  as  the  larger  of  the  two. 
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III.  MULTIPLEXING  IN  THE  IMPLEMENTATION  OF  DIGITAL  FILTERS 


In  a sense,  any  time  a computer  containing  a single  multiplier  and 
adder  is  used  to  implement  a digital  filter  multiplexing  is  done.  To  see 
this  simply  consider  the  nonrecursive  digital  filter  given  by 
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a,  X 

k n-k 


(72) 


If  a single  multiplier  is  employed  each  product  in  (72)  is  computed 
sequentially  and  multiplexing  of  the  multiplier  is  effectively  done  in 
time.  Wlien  more  than  one  filter  is  Implemented  on  a digital  computer 
time  division  multiplexing  can  be  employed.  This  is  shown  in  Fig.  2 
where  two  channels  are  handled  by  a single  computer.  At  time  t^Q  a 
sample  is  taken  on  channel  1 and  an  output  generated  at  time 
computational  time  is  less  than  T/2  seconds.  At  time  t^Q  a sample  is 
taken  from  channel  2 and  the  output  is  computed  for  use  at  time  *^21 ' 

If,  however,  there  is  ample  computational  time,  both  samples  may 
be  taken  at  time  t^  . The  response  to  each  input  sample  is  then  com- 
puted sequentially  and  both  outputs  are  generated  at  time  t^  for  use 
by  the  separate  channels.  This  is  shown  in  Fig.  3.  Here,  at  the  samp- 
ling Instants  each  input  is  sampled  by  a sample-and-hold  circuit  and 
the  inputs  are  converted  into  binary  and  stored  in  memory.  The  output 
for  the  first  filter  is  computed  during  Computer  Time  1 and  the  output 
for  the  second  filter  is  generated  during  Computer  Time  2.  The  two 
outputs  are  stored  in  buffer  registers  that  are  clocked  to  the  output 
channel  lines  at  the  end  of  tlie  sampling  interval  . 

At  this  point  it  is  important  to  realize  that  most  computer  .md/or 
data  acquisition  system  manufacturers  do  not  design  the  control  for  the 
analog-to-digital  conversion  and  digital-to-analog  conversion  subsystems 
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in  the  right  way.  Since  the  coefficients  and  thus  the  digital  filter 
characteristics  are  very  dependent  on  the  sampling  interval  T the 
control  data  acquisition  system  must  be  carefully  designed.  Thus,  the 
real  time  clock  should  be  an  Integral  part  of  the  data  acquisition  system 
and  control  should  not  be  done  via  interrupts. 

For  example,  in  many  computer  systems  the  real  time  clock  is  a 
feature  of  the  computer.  This  consists  of  a preset  counter  and  when 
the  count  reaches  zero  an  interrupt  is  activated  to  tell  the  computer  to 
take  a sample  from  the  data.  This  is  done  via  a subroutine  that  controls 
the  analog-to-digital  converter.  However,  in  normal  operation,  the  time 
to  jump  to  the  service  routine  may  vary  and  in  many  cases  the  variation 
will  be  a large  percentage  of  the  sampling  interval.  As  a result  the 
actual  sampling  interval  is  varying  in  a random  fashion  or  is  in  error 
by  a constant.  The  characteristics  of  the  filter  will  now  deviate  from 
the  desired  characteristics.  For  example,  if  a band-stop  filter  is  being 
implemented,  the  center  frequency  will  be  shifted.  It  is  also  possible 
to  achieve  instability  in  the  digital  filter. 

In  a good  system,  a programmable  real-time  clock  will  be  specified 
as  part  of  the  data  acquisition  system.  This  clock  will  clock  a sample 
point  from  the  analog-to-digltal  converter  and  then  interrupt  the  computer 
to  indicate  a sample  has  been  taken  and  is  ready  for  processing.  If  a 
digital- to-analog  converter  is  used,  the  filter  output  will  be  computed 
and  stored  in  a buffer.  At  the  time  when  a new  sample  is  taken,  the 
output  data  will  be  transferred  to  the  digital-to-analog  converter. 

Wfien  more  than  one  channel  is  serviced  by  the  computer,  a multi- 
plexer is  used  in  front  of  the  analog-to-dlgital  converter.  Here  care 
must  be  taken  to  avoid  a time  lag  between  samples.  This  can  be  handled 


by  using  a samp] e-and-hold  on  each  channel  before  the  multiplexer. 

Another  alternative  is  to  allow  a fixed  delay  between  sampling  points. 

By  including  tlie  same  delay  at  the  output  the  sampling  interval  for  a 
given  channel  will  be  constant. 

If  a computer  is  used  to  implement  the  digital  filters  in  an 
aircraft  control  and  avionics  system,  the  data  will  be  transmitted  to  and 
from  the  computer  in  two  ways.  First  for  every  channel,  an  indiyidual 
wire  can  be  used.  Secondly,  a bus  structure  can  be  employed.  In  this 
case,  the  data  will  be  time  multiplexed  on  the  bus.  If  a single  bus 
is  employed,  both  the  input  and  output  data  will  use  the  bus. 

lilhen  an  individual  wire  is  used  for  each  channel,  the  synchronizing 
of  the  data  to  achieve  a uniform  sampling  time  is  not  difficult  and  the 
data  acquisition  system  should  be  specified  as  outlined  previously, 
llhen  a bus  structure  is  used,  unless  the  bus  is  operated  synchronously 
there  is  a danger  of  introducing  variations  in  the  sampling  interval 
for  a given  channel.  These  variations  can  cause  the  digital  filter 
to  perform  improperly.  As  a result,  the  designer  must  carefully  specify 
the  amount  of  allowable  error  in  the  sampling  interval  to  establish  a bus 
operating  procedure  that  will  satisfy  the  sampling  Interval  specification. 

In  terms  of  reliability,  redundancy  is  usually  employed  by  using 
both  redundant  computers  and  redundant  bus  structures.  In  actual  opera- 
tion faults  are  detected  using  a majority  vote. 

As  the  Integrated  circuit  manufacturers  develop  the  new  microprocessor 
systems,  a more  realistic  procedure  for  handling  the  digital  filtering 
operations  consists  of  a distributed  computing  system.  In  such  a system 

i 

a separate  microprocessor  can  be  employed  to  implement  each  digital 
filter.  By  placing  each  microprocessor  under  the  control  of  the  central 
computer  rapid  reprogramming  of  the  filter  coefficients  can  be  achieved 
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as  operating  conditions  change.  Improved  reliability  in  terms  of 
aircraft  operation  can  also  be  achieved  by  a filter  priority  structure. 
This  is  illustrated  in  Fig.  4.  Referring  to  the  figure,  three  micro- 
processors are  used  to  implement  three  digital  filters  in  channels  1, 

2 and  3.  The  priority  has  been  predetermined  to  be  channel  1 highest 
priority,  channel  2 next  highest  priority  and  channel  3 lowest  priority. 
The  channel  inputs  and  outputs  are  switched  into  the  appropriate  micro- 
processors under  control  of  tlie  central  computer.  If  microprocessor  1 
fails,  channel  1 is  switched  into  microprocessor  3 and  channel  3 is  dis- 
connected. Microprocessor  3 is  then  reprogrammed  to  handle  the  filter 
for  channel  1.  If  desired,  additional  redundancy  can  be  Included  in 
each  microprocessor. 

If  a single  data  bus  is  employed,  the  control  computer  will  program 
the  proper  address  into  each  microprocessor  so  that  the  data  will  be 
pointed  into  the  correct  microprocessor.  In  the  event  of  a high  priority 
microprocessor  failure,  the  control  computer  would  reprogram  the 
address  into  a low  priority  microprocessor  and  discontinue  pointing  data 
from  the  low  priority  channel. 

Because  the  present  structure  of  a digital  flight  control  and 
digital  avionics  systems  appears  to  be  based  on  using  one  computer  for 
handling  all  of  the  digital  filtering  operations,  the  development  of 
simulation  software  for  multiple  digital  filters  for  digital  control 
was  started  as  a part  of  the  project.  The  computer  was  a Digital 
Equipment  Corporation  PDP-11.  The  complete  description  of  this  work  was 
sent  to  the  contract  monitor  as  a thesis  referenced  in  [14].  In  this 
report  only  the  necessary  details  are  included. 

The  general  structure  of  the  multiplexed  digital  control  system 
is  shown  in  Fig.  5 for  four  individual  control  channels.  In  tlie 
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simulation,  each  channel  input  is  sequentially  multiplexed.  When  a 
sample  is  taken  from  channel  1,  the  proper  filter  algorithm  is  called 
and  the  filter  output  is  computed  and  sent  to  the  digital-to-analog 
converter.  Then  channel  2 is  sampled  etc.  Thus,  each  channel  is 
handled  separately  with  the  proper  filter  coefficients  being  placed  in 
a memory  stack  corresponding  to  the  given  channel.  The  flow  chart  is 
shown  in  Fig.  6. 

The  software  package  for  this  simulation  is  not  complete  because 
of  the  change  in  the  data  acquisition  system.  However,  the  program 
using  an  outdated  data  acquisition  system  is  complete  and  is  described 
in  Reference  [14]. 

To  summarize,  the  current  trend  in  the  digital  avionics  and  flight 
control  system  is  to  utilize  redundant  computers.  This  implies  that  a 
single  computer  will  be  used  to  handle  all  of  the  digital  filtering 
operations  in  the  aircraft.  The  data  will  be  handled  via  a bus  and 
it  will  be  necessary  to  design  this  in  such  a way  that  the  sampling  of 
data  will  correspond  to  the  way  that  sampling  is  done  in  the  filter 
design  process. 

The  multiplexing  for  this  configuration  will  consist  of  a single 
computer  used  to  Implement  a number  of  filters.  For  a given  number  of 
filters,  all  of  the  input  sample  points  at  a given  time  can  be  taken  at 
once  and  all  filter  output  computed  during  a single  sampling 
interval  or  a single  filter  can  be  implemented  during  a sampling 
interval.  In  either  case  the  data  acquisition  system  must  be  carefully 
specified  and  designed.  If  this  is  not  done  timing  will  be  difficult 
and  a considerable  software  overhead  will  be  required  to  acquire  and 
send  out  data. 
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In  future  designs  it  is  desirable  to  consider  the  use  of  distributed 
computing  systems  with  microprocessors  being  used  to  implement  digital 
filtering.  By  employing  a priority  multiplexed  switching  system  under 
control  of  a central  computer,  the  highest  priority  filtering  operations 
can  be  maintained  during  flight. 
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IV . CONCLUSIONS  AND  RKCOMMENUATIONS  FOR  FUTURE  WORK 


Currently  it  appears  the  digital  filter  design  using  the  bilinear  z 
transform  is  most  suitable  for  digital  flight  control.  However  for  many 
avionics  applications,  time  domain  synthesis  may  be  more  effective.  In 
the  future  it  is  important  that  attention  be  given  to  developing  design 
methods  that  are  more  suited  to  the  z domain.  For  example,  the  sampling 
interval  should  be  used  as  a design  parameter.  For  effective  design  it 
is  necessary  to  develop  good  interactive  software  to  work  with  a graphics 
terminal . 

In  addition  it  is  Important  that  the  designer  carefully  consider  the 
implementation  of  a digital  filter  when  doing  the  design.  Specifically 
the  effects  of  rounding  must  be  considered.  By  the  use  of  interactive 
design  software,  and  the  procedures  developed  in  Appendices  3,  4,  and  5, 
the  word  length  needed  for  an  implementation  can  easily  be  determined. 

When  the  design  for  the  complete  aircraft  system  is  developed, 
the  method  of  moving  data  within  the  aircraft  must  be  carefully  considered. 
If  this  is  not  done  the  shift  in  sampling  interval  can  cause  severe  prob- 
lems in  the  way  a digital  filter  performs.  As  part  of  this  work, 
attention  must  be  given  to  the  specification  and  design  of  the  data 
acquisition  systems  in  order  to  minimize  software  overhead  and  to 
insure  proper  filter  operation. 

While  multiplexing  can  be  employed  future  systems  should 
utilize  a distributed  computing  structure  for  greater  reliability.  This 
will  reduce  tlie  need  for  a large  central  computer  and  will  allow  greater 
flexibility  in  the  implementation  of  digital  filtering  algorithms. 

Last,  the  need  for  developing  good  interactive  graphics  is  again 
emphasized.  In  addition,  a simulation  system  should  be  establ  islied  to 

34 


i 

I 


allow  Air  Force  personnel  to  carefully  develop  specifications  for  digital 
avionics  iind  flight  control  systems.  The  effective  use  of  simulation  and 
interactive  design  will  result  in  more  effective,  less  expensive  and 
more  reliable  aircraft  systems. 
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Fi".  3 Block  Diagram  for  a Two  Channel  System 
lirtiere  Both  Kilter  Oiiti)uts  are  Cenerated 
in  One  Sampling  Interval 
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INTRODUCTION 


This  report  contains  the  documentation  for  the  BRASS  program.  It 
consists  of  the  design  procedure  used,  a description  of  the  program, 
and  design  examples  using  the  program. 

The  purpose  of  the  BRASS  program  is  the  design  of  either  a 
maximally  flat  Butterworth  or  a Jhebychev  filter  with  equal  ripple  in 
the  pass  band.  For  each  type  of  filter  there  is  a choice  of  band-pass 
or  band-stop  filters.  Starting  with  an  analog  filter,  the  bilinear 
Z transform  is  used  to  design  an  equivalent  digital  filter.  The  user 
enters  the  low-pass  filter  order,  the  type  of  filter  desired,  the 
sampling  interval,  the  upper  and  lower  cutoff  frequencies,  the 
starting  frequency  and  frequency  increment,  and  if  a Chebychev  filter 
is  being  designed,  the  ripple.  The  low-pass  filter  sections  are 
transformed  to  second  order  band-pass  or  band-stop  sections.  Then  the 
program  generates  the  digital  filter  coefficients  for  up  to  six 
second  order  sections  in  cascade  or  up  to  a 12th  order  filter.  The 
design  is  carried  out  in  the  frequency  domain.  The  program  calculates 
the  transfer  function  coefficients  for  each  second  order  section,  the 
magnitude  function  for  each  section,  and  the  final  cascaded  filter 
magnitude  response  over  the  frequency  interval  specified  by  the  input. 

The  BRASS  program,  written  in  Fortran  IV  is  supplied  as  a card 
deck  with  this  report.  The  program  is  in  the  form  of  a subroutine  and 
can  be  used  as  is  by  a call  statement  from  the  main  program.  Data 
may  be  input  via  cards  with  output  available  through  a line  printer. 
The  input/output  devices  may  be  altered  as  explained  in  this  report. 
Graphic  routines  may  easily  be  appended  to  the  program. 


» 

'i 

t 
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Design  Procedure 


I . 

A.  Preliminary  Discussion 

One  common  method  of  designing  a digital  filter  is  to  start 
with  an  anaiug  transfer  function  H(S)  and  transform  it  to  the 
digital  transfer  function  H(Z). 

The  transfer  function  of  a second  order  digital  filter  in  the 
Z domain  is  given  by 

K. (A„Z^  + A.Z  + A^) 

H(Z)  = i ^ (1) 

Z + B^Z  + B, 

where  the  A's  and  B's  are  the  coefficients  of  the  numerator  and 
denominator  respectively.  This  program  will  calculate  the  scale 
factor  and  the  coefficients  A^,  A^ , A^ , B^^,  and  B^ . The 

transformation  used  is  the  extended  bilinear  Z transform 


S 


2 Z - 1, 

t4  + 1' 


(2) 


where  T is  the  sampling  interval.  When  the  extended  bilinear  Z 
transform  is  employed,  the  desired  frequencies  must  first  be  pre- 
warped to  make  them  compatible  with  the  digital  filter.  In  the 

band-pass  and  band-stop  filters,  the  upper  and  lower  cutoff  frequencies 

* 

and  the  center  frequency  of  the  filter  are  of  Interest.  Calling  the 
upper  and  lower  frequencies  u)^  and  respectively',  the  pre- 

warped upper  (WDU)  , lower  (i'lDL)  , and  center  (WDM)  frequencies  and  the 
bandwidtli  between  WDU  and  WDI  are  found  by 


ry  U»  X 

WDU  = ^ tan(-^) 

WDL  = I tan(-^) 

2 

» WDM  = Y tan 

WB  = WDU  - WDL 

oj^  and  are  specified  by  the  designer  and  the  prewarping  is  done 

by  the  program. 

In  the  design  procedure  for  all  band-pass  and  band-stop  filters  of 
order  n',  (n'  even),  the  program  begins  by  first  finding  the  poles 

for  the  corresponding  n'/2  order  low-pass  filter.  The  low-pass 
filter  is  then  transformed  into  a band-pass  or  band-stop  filter  of 
order  n'.  (n'  = 2n) . 

B.  Butterworth  Band-Pass  Filter 

We  start  with  a normalized  second  order  low-pass  Butterworth 
filter  transfer  function  in  the  S plane 

H(S)  = (4) 

S + 2Scos6  + 1 

where  the  angle  9 is  in  degrees  (in  the  program)  and  may  be  found 
fi am  the  Butterworth  circle  and  the  relationship 

S = (5) 

where  n is  the  order  of  the  low-pass  filter  and  m=  1,  2,...,  n. 

This  relationship  is  determined  by  the  following  procedure.  By 
definition,  a filter  is  nth  order  Butterworth  low-pass  if  its  gain 
characteristic  is 
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1 


Ill^Cja.) 


1 ^ 

c 


2n 


(6) 


wliere  a is  the  DC  gain,  u)^  is  the  desired  cutoff  frequency  and  n 
is  the  order  of  the  low-pass  filter. 

In  the  design,  the  poles  of  H(S)  must  be  found.  The  procedure 
is  as  follows: 


|H^(jw)|  = H^(jw)H^(ju)  = ( joj)H^  (-jcij) 


= [H(S)H(-S)] 


S = jw 


1 + (— ) 

u 

c 


2n 


2n 


S 

lo)  = — 
1 


1 + (t^) 

J“c 


1 + 


1 + 


1 - 


-2 


, for  n even 


, for  n odd 


(7) 


Setting  the  denominators  equal  to  zero. 


Thus,  the  pole  locations  are  the  2n  roots  of  ±1,  depending  on 
whether  the  low-pass  filter  order  Is  odd  or  even.  These  roots  are 
located  on  a circle  with  radius  w centered  at  the  origin  of  the  S 


i 


U) 

c 


(8) 


^1 


plane  and  have  symmetry  with  respect  to  both  real  and  imaginary  axes. 
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For  n odd,  a pair  of  roots  are  on  the  real  axis  and  the  rest  are 


separated  by  n/n  radians.  For  n even,  a pair  of  roots  are  located 
Tr/2n  radians  from  the  real  axis  and  the  rest  are  again  separated  by 
Ti/n  radians.  No  roots  are  on  the  imaginary  axis,  for  either  even  or 
odd  n . 


Let  Pj^,...,p^^  be  the  roots.  From  the  symmetry  of  the  pole 
locations,  if  p^^,...,p^^  are  the  roots  lying  in  the  right-half  plane, 
the  left-half  plane  roots  are  -p^,...,-p^.  The  magnitude-squared 
function  can  then  be  written  as 


2.  ,n  2n 
a (-1) 

H^(S)Hn(-S)  = (s  + p^...(s  + p_^)(S  - p,  ) . . . (S  - p„) 


(9) 


n 


I' 


n 


To  be  stable,  H (S)  must  have  all  its  poles  in  the  left-hand  plane, 
n 

thus 


SLb) 

" Ts  + p.). . .(s  + p„) 


(10) 


The  program  is  written  with  unity  gain  at  DC,  (cu  = 0) , therefore 
a = 1. 

In  order  to  locate  the  poles  as  specified  above,  consider  the 
following  sot  of  equations. 


1 = 

^ijir(2m  - 1) 

, m = 1 , 2 , . . 

. , n ; for  n even 

1 = 

_^ij2TTk 

o 

II 

, . , n ; for  n odd 

(11) 


Substituting  equations  (11)  into  equations  (8)  yields 


S , ijTT(2m  - l)/2n  , „ 

— =_gj'  m=l,2,..,n;  for  n even 


0)  -tm 
c 


(12) 


= - 


+ juk/n 


, k = 0,  1,...,  n;  for  n odd 


Equations  (12)  will  give  the  pole  locations  as  described  above. 
Consider  the  form  of  equations  (12) 

S = -0)^  v ■ -J  = j^[-cos9  ♦ jsin9]  . (13) 


From  this  relationship,  it  can  been  seen  that  the  magnitude  for  each 

pole  is  , regardless  of  the  angle,  and  thus  all  the  poles  lie  on  a 

circle  with  radius  w . 

c 

As  an  example,  consider  a second  order  filter,  n = 2. 


[.§_]  = _gij7T(2m  - l)/4 

w +m 
c 


m = 1 , 2 


S 


±1 


0)  /±45“ 
c 


S,,  = u /±135° 

±1  c 

6 = 45° 

The  relationship  of  these  roots  about  the  circle  of  radius  is 

illustrated  in  Figure  1.  The  angle  6 is  always  measured  from  the 
negative  real  axis. 

In  the  program,  only  the  angle (s)  less  than  90°  are  considered 
so  that  the  poles  lie  in  the  left-half  plane  because  poles  in  the 
left-half  plane  are  stable.  Putting  0 = 45“  into  equation  (4) 
yields  poles  at  -0.707  ±j0.707.  These  locations  are  in  the  left-half 
plane.  From  equations  (12),  for  low-pass  filter  orders  n = 1,  2,..., 
6,  the  values  of  9 are  given  below. 


i 
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Second 

Band-Pass 

■ow-Pass 

Order 

Band-Stop 

Filter 

Cascaded 

Filter 

Order 

Angle 

Sections 

Order 

n 

0 

N 

n' 

1 

0° 

1 

2 

2 

45“ 

2 

4 

3 

60“ , 0“ 

3 

6 

4 

22.5“  , 67.5“ 

4 

8 

5 

72“  , 36“  , 0“ 

5 

10 

6 

75“  , 45“ , 15“ 

6 

12 

n is  the  order  of  the  low-pass  filter  and  is  used  to  determine  pole 
locations,  n is  also  the  number  of  second  order  band-pass  or  band- 
stop  sections  which  results  from  the  transformation  of  the  low-pass 
filter  sections  and  which  will  be  cascaded  to  form  the  band-pass  or 
band-stop  filters  of  order  n'.  The  transformation  is  explained  below. 
The  calculated  angles  are  incorporated  in  the  program  in  the  order 
given  above. 

Given  the  normalized  second  order  low-pass  transfer  function 
equation  (4),  we  transform  this  low-pass  into  a band-pass  transfer 
function  for  some  bandwidth  WB , and  center  frequency  WDM  by 
using  the  transform 


2 2 
^ S + WDM 

^ " SWB 


(14) 


Equation  (4)  then  transforms  to  a 4th  order  transfer  function 
H(S)  = 


2 2 
S WB 


s'^  + 2WBCOS0  + S^(2WDM^  + WB^)  + S2WB  WDM^  cos0  + WDM^ 


(15) 


Using  the  root  finding  subroutine  "POLRT"  from  tlie  IBM  Scientific 
Subroutine  Package  (SSP)  , the  roots  of  the  denominator  of  equation  (15) 
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are  found.  (Note:  POLRT  has  been  attached  to  BPASS  as  a double 


precision  subroutine  and  is  included  in  tlie  card  deck).  The  roots 
found  will  be  complex  conjugate  pairs.  Calling  the  real  and  imaginary 
parts  of  the  pairs  RE^ , AIM^ , RE^,  AIM^  equation  (15)  is  factored 
to  yield  two  cascaded  second  order  sections 


H(S)  = 


SWB 


SWB 


2SRE^  + RE^  + AIM^ 


- 2SRE2  + RE^  + AIM^ 


(16) 


For  each  0 of  a given  N,  the  program  calculates  roots  for  both 
sections  of  equation  (16)  and  labels  them  the  ith  and  the  ith  + 1 
section.  If  N,  the  number  of  second  order  sections  specified,  is 
even,  the  program  will  calculate  N pairs  of  RE  and  AIM  values 
or  2N  = n’  roots.  If  N is  odd,  the  last  value  of  9 is  0. 
Substituting  0=0  into  equation  (4)  and  factoring  yields  two 
identical  first  order  sections,  1/(S  + 1).  The  program  will  calculate 
N + ' pairs  of  RE  and  AIM  values,  but  because  the  last  two  pairs 
are  the  same  due  to  the  identical  first  order  sections,  the  last  pair 
will  not  be  used. 

Because  both  second  order  sections  of  equation  (16)  are  of  the 
same  format,  we  will  deal  with  only  one  section,  the  ith  section 
and  let 


-2RE.  = D. 

^ ‘ ( 

RE^  + AIM^  = C. 

1 11 

The  design  of  an  n'th  order  band-pass  or  band-stop  filter  leads  to 
n'/2  second  order  sections.  Substituting  equations  (17)  into  (uie 
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section  of  equation  (16)  yields  the  transfer  function  for  the  ith 
section 


H.  (S) 


SWB 

S“  + SD.  + C. 


(18) 


1 1 

The  extended  bilinear  Z transform,  equation  (2)  is  used  to  get  to 
the  digital  domain.  Employing  equation  (2)  on  equation  (18)  yields 
li^(Z)  for  the  ith  second  order  section. 


H.  (Z) 


2 2 2 
^ 


9 / 2D. 

Z^  (—  + ^ 

^^2  T 


+ C ) + Z(2C  - 

T T 


2D  . 
1 


(19) 


+ C.) 


Putting  the  denominator  of  equation  (19)  in  monic  form  yields  the 
transfer  function  for  the  1th  second  order  stage  of  the  filter 


H.(Z) 


K^.(AqZ^  + A^Z  + A^) 

Z^  + B, .Z  + B„. 

li  2i 


(20) 


This  equation  is  the  same  as  equation  (1)  with  the  exception  of  the 
subscripts.  For  all  four  filter  types  discussed  here,  the  scale 
factor,  Kj^,  and  coefficients  B^  and  B2  are  a function  of  the 
section  calculated,  wliile  the  coefficients  Aq,  A^,  and  A2  are  the 
same  for  all  sections  calculated.  In  going  from  equation  (19)  to 
equation  (20)  we  have 
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Aq  = 


= 0 


A2  = - 

, 2D. 

1 


(21) 


2C 


11 


i_i! 

G. 

1 


/ 2D. 

-j  - + C. 

t2  T 


Letting  Z = e - for  S = jio  and  taking  the  magnitude  of 

H . ( jw)  we  have 


|H. (jw) I = 


^AqCOs(2uT)  + A cos(wT)  + A„)^  + (A„sin(2ajT)  + A,  sin(tuT))^ 


/(cos(2o)T)  + B^^cos(wT)  + B.,^)  + (sln(2uT)  + B^^sin(u)T) )‘ 

(22) 

The  magnitude  function,  equation  (22)  is  the  same  for  all  the  filters 
discussed  in  this  report. 

C.  Butterworth  Band-Stop  Filter 

The  design  procedure  is  almost  exactly  the  same  as  that 
of  the  Butterworth  band-pass  filter,  except  that  the  transformation 
to  band-stop  is  the  reciprocal  of  equation  (lA),  i.e. 


S ->■ 


SWB 


2 2 
S + WDM 


(23) 
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(24) 


and  we  find  H^(S)  to  be 


H,(S) 


2 2 
S + WDM 


S + SD.  + C. 

1 1 


After  employing  the  extended  bilinear  Z transform,  equation  (2), 
we  have 


A = = ^ + WDM^ 

T 

A = 2WDM^  - ~ 
i 


(25) 


and  B , B„ . , K are  the  same  functions  of  C.  and  D.  as  in 

iizili  1 1 

equation  (21).  These  coefficients  are  then  used  in  the  calculation 
of  equation  (22)  to  find  |n.(ju))|. 

D.  Chebychev  Band-Pass  Filter 

The  Chebychev  filter  ripples  with  equal  amplitude  in  the 
pass-band.  The  amount  of  ripple  is  specified  by  the  quantity  6 
(labeled  RIP  in  the  program).  The  poles  of  the  filter  are  found  on 
an  ellipse  described  by  two  Butterworth  circles  of  radii  A and  B 
with  A < B.  The  location  of  tlie  poles  on  the  ellipse  is  a function  of 
the  ripple  and  is  given  by  the  following  equation: 


B.  A = t((vi  ^ + 1 + t 


« . (,f  2 , 1 „ (26) 


where 


(27) 


and  N is  numerically  equal  to  the  order  of  the  low-pass  filter 
which  is  transformed  to  yield  the  band-pass  filter.  B is  given 
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for  the  plus  sign  and  A for  the  minus  sign.  The  Chebychev  ellipse 
then  has  major  axis  B and  minor  axis  A.  The  location  of  the  S 
plane  poles  on  the  ellipse  is  given  by 


Real  Part  = A cos0 
Imaginary  Part  = sin6 


(28) 


The  9's  are  the  same  as  given  for  the  corresponding  order 
Butterworth  filter.  An  example  of  Chebychev  pole  locations  is 
illustrated  in  Figure  2.  For  ^ B = 1 in  a fourth  order 

filter,  9 = 22.5°  and  67.5°.  The  Chebychev  pole  locations  are 
determined  from  equations  (26),  (27)  and  (28). 

The  analog  second  order  Chebychev  low-pass  filter  is 


H(S)  = 


^2 

1 

2/N 

U1  + e^J 

(29) 


S + KgS  K, 


e is  calculated  from  equation  (27)  and  N is  equal  to  the  order  of 
the  low-pass  filter  which  is  transformed  to  yield  the  band-pass  filter 
Kg  and  K^  are  calculated  by 


K_  = 2Acos9 
o 

2 2 2 2 
K2  = A cos  9 + B sin  9 


(30) 

(31) 


The  substitution  of  the  low-pass  to  band-pass  transformation, 
equation  (14),  into  equation  (29)  yields 


H(S)  = - 


2 2 
S WB 

1 

2/N 

,/l  + e^. 

+ S^K„WB  + S^(2WDM^  + K„WB^)  + SK.WDM^WB  + WDM'^ 

O 4 0 


(32) 
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After  finding  the  roots  of  equation  #32)  and  making  the  substitutions 
given  by  equations  (17)  we  find  the  ith  second  order  section 

SWK 

H.  (S)  = ^ 

^ S“  + SD.  + C. 

1 L 


K3  = 


' i + c 


1/N 


(33) 


Applying  the  extended  bilinear  Z transform  equation  (2)  yields  an 
equation  of  the  form  of  equation  (20)  where 


B, . and  B_ . are  the  same  functions  of  C.  and  D.  given  by 
li  2i  11'’^ 

equations  (21).  These  coefficients  are  then  used  in  equation  (22)  to 
f ind  1 H3 ( jw) I . 

E.  Chebychev  Band-Stop  Filter 

Given  equation  (29)  for  H(S)  we  apply  the  low-pass  to 
band-stop  transformation  equation  (23)  to  obtain  the  4th  order 
transfer  function 

2 2 2 
(S  + WDM  ) K2 

H (S)  = 3 2 2 2 2 4 ' 

K-S  + S K-WB  + S (WB  + 2K„WDM  ) + SKQl>mM  WB  + K„WDM 

Z o Z o 4- 

N is  equal  to  the  order  of  the  low-pass  filter  which  is  transformed 
to  yield  the  band-stop  filter. 


U' 


+ c J 


2/N 


1 


I 

.» 

j 

r 

, 
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After  finding  the  roots  of  equation  (35)  and  making  the 


substitutions  given  by  equations  (17)  tiie  ith  second  order  section  is 


(S^  + WDM")K. 


Hi(S)  = 


S + SD  + C. 

1 1 


S • 'S 


/I  + e ' 


i/N 


(36) 


Applying  the  extended  bilinear  Z transform  equation  (2)  yields  an 
equation  of  the  form  of  equation  (20)  where 


*0  ■ *2  ■ ^ 


A = 2WDM^  - ^ 
1 


(37) 


11  G. 
1 


B, . and  B„ . are  the  same  functions  of  C.  and  D.  given  by 
li  2i  11 

equations  (21).  These  coefficients  are  then  used  in  equation  (22) 
to  find  I H.  (juj)  | . 
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II . Using  the  ProRram 

The  first  data  card  read  into  the  program  contains  the  number  of 
second  order  sections  to  be  cascaded,  N,  and  tlie  type  of  filter 
desired,  KN.  N is  equal  to  1,  2,...,  or  6,  which  corresponds  to 
the  order  of  the  low-pass  filter,  and  hence  corresponds  to  the  2nd, 
4th,,..,  or  12th  order  band  pass  or  band  stop  filter  respectively. 

KN  is  the  type  of  filter  desired.  The  values  of  KN  specifies 
one  of  the  four  choices  given  by 


KN 


Type 


Butterworth  Band-Pass 
Butterworth  Band-Stop 
Chebychev  Band-Pass 
Chebychev  Band-Stop 


The  format  on  the  N,  KN  card  is  212. 

The  second  data  card  read  in  is  the  sampling  interval  T in 

F10.6  format.  When  choosing  T,  l/T  should  be  approximately  equal 

to  ten  times  the  center  frequency  (WDM). 

The  tliird  data  card  read  in  contains  the  values  of  the  upper 

and  lower  cutoff  frequencies,  m and  w,  , in  2F10.4  format.  For 

u 1 

the  Butterwortli  filters,  the  cutoff  frequencies  are  the  -3db  cutoff 
frequencies.  For  the  Chebychev  filters,  the  magnitude  of  the 
response  i s I / ( 1 -t-  ) “ I - 6 at  the  cutoff  frequencies,  oj  is  in 

radians.  is  t lie  ripple  factor. 

If  the  desired  filter  is  Chebychev,  i.e.,  KN  = 3 or  4,  the 
next  data  card  caintains  the  ripple  (RIP)  factor  in  F5.3  format. 

If  the  desired  filter  is  Butterworth,  i.e.,  KN  = 1 or  2,  this 
card  is  omitted  from  the  data  deck. 


I 
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The  final  data  card  is  the  starting  frequency  (FREQl)  and  the 
frequency  increments  (DELT)  in  radians.  The  format  of  the  FREQl, 
DELT  card  is  2F10.4.  Determine  DELT  by  the  following: 


DELT  = 


final  frequency  - starting  frequency 


This  is  necessary  because  there  are  1024  frequency  data  points 
calculated  in  the  program.  Choose  FREQl  and  DELT  to  insure  that 
calculated  values  will  include  the  data  of  interest.  For  maximum 
efficiency  of  the  program,  DELT  should  be  a multiple  of  2 so  no 
decimal  to  binary  conversion  errors  are  incurred. 

The  digital  filter  coefficients  are  computed  and  printed  out  for 
each  second  order  section.  The  full  filter  magnitude  response,  as  well 
as  each  section  magnitude  response,  is  printed  for  each  of  the 
specified  frequency  increments.  When  there  is  only  one  second  order 
section,  the  section  magnitude  response  is  the  full  filter  magnitude 
response  and  is  only  printed  once. 

The  program  may  be  easily  modified  to  incorporate  a graphics 
display  of  the  magnitude  response.  There  is  a comment  card  in  the 
BPASS  program  indicating  where  the  graphics  subroutine  call  card 
should  be  inserted. 

The  program  is  written  with  input  obtained  via  device  4 and 
output  written  to  device  6.  These  numbers  should  be  assigned  to  the 
appropriate  devices  prior  to  running  the  program. 

The  program  was  developed  on  a PDP-11/20  with  a DOS/BATCH 
operating  system.  Trial  runs  frequently  used  a TTY  terminal  as  well 
as  a card  reader  for  input  (device  4);  and  a TTY  terminal  as  well 


as  a line  printer  for  output  (device  6).  Double  precision  aritlimetic 
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is  employed.  To  decrease  required  memory  storage,  only  the  frequency 
interval  values  and  the  full  magnitude  response  are  saved.  The 
section  magnitude  responses  are  printed  out,  but  are  not  stored. 

The  program  will  produce  approximately  21  pages  of  output. 

Shown  below  are  sample  deck  set-ups. 


Format 


Examp 1 e 


212 

F10.6 

2F10.4 

F5.3 

2F10.4 


0504  (5  sections  Chebychev  band-stop) 

0.002  (T  = 0.002) 

60  40  (u)  = 60,  lii,  = 40  radians) 

u 1 

0.10  (Ripple  amplitude  = 0.10) 

0 0.1  (Start  at  w = 0.  Steps  of 
0.1  radian.  Will  finish 
just  past  (j  = 102  radians.) 


212 

F10.6 

2F10.4 

2F10.4 


0401  (4  sections  Butterworth  band-pass) 
0.002  (T  = 0.002) 

60  40  (u)  = 60,  = 40  radians) 

0 0.1  (Start  at  ij  = 0.  Steps  of 
0.1  radian.  Will  finish 
just  past  ui  = 102  radians)  . 


The  following  pages  contain  annotated  examples  of  output  data. 
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This  is  an  example  of  the  output  for  an  8th  order  Butterworth 

band-stop  filter  (0402)  with  T = 0.002,  ui  = 60  radians,  and 

u 

01^  = 40  radians.  The  starting  frequency  is  0 radian  and  the  frequency 
increment  is  0.1  radian. 


WDU  = 60.07210  WDL  = 40.02135  WDM  = 49.02902  WB  = 20.05076 
T = 0.20000E-02 


THE  ROOTS  OF  THE  FILTER  ARE  GIVEN  BELOW 
REAL(l)  = -8.52659476  IMAGINARY(l)  = -44.46786740 
REAL(2)  = -9.99788916  IMAGINARY(2)  = -52.14095930 
REAL(3)  = -4.55076557  1MAGINARY(3)  = -59.01588870 
REAL(4)  = -3.12232691  IMAGINARY(4)  = -40.49140500 


THE  GOEFFICIENTS  OF  EACH  DIGITAL  FILTER  SECOND  ORDER  SECTION  ARE 
GIVEN  BELOW 


FOR 

I = 1 

= 

0.10024038E+07 

= 

-0.19951923E+07 

= 

0.10024038E+07 

= 

0.98125481E-06 

= 

-0.19584863E+01 

®2 

= 

0.96653295E+00 

FOR 

1 = 2 

A° 

= 

0.10024038E+07 

^1 

-0.19951923E+07 

= 

0.10024038E+07 

= 

0.97769448E-06 

= 

-0.19498774E+01 

^2 

= 

0.96090048E+00 

FOR 

I = 3 

0.10024038E+07 

- 

-0.19951923E+07 

= 

0.10024038E+07 

•"I 

0.98755180E-06 

= 

-0.19681836E+01 

^2 

= 

0.98202353E+00 

FOR 

II 

'2 

0.10024038E+07 

-0.19951923E+07 

= 

0.10024038E+07 

®2 

= 

0.99216788E-06 

^1 

= 

-0. 198106 30E+01 

= 

0.98760851E+00 

W H 

0.0000  O.lOOOOE+01 
0.1000  O.lOOOOE+01 
0.2000  O.lOOOOE+01 
0.3000  0.99999E+00 

0.4000  O.lOOOOE+01 
0.5000  O.lOOOOE+01 


HI 

0.11726E+01 

0,11726E+01 

0.11726E+01 

0.11726E+01 

0.11726E+01 

0.11726E+01 


H2 

0.85284E+00 
0.852 84 E+00 
0. 85284 E+00 
0.85283E+00 
0.85283E+00 
0.85282E+00 


H3 

0.68611E+00 

0.68611E+00 

0.68611E+00 

0.68610E+00 

0.68609E+00 

0.68609E+00 


H4 

0.14575E+01 

0.14575E+01 

0.14575E+01 

0.14575E+01 

0.14576E+01 

0.14576E+01 
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r : r 

I 

WDU  is  the  prewarped  upper  frequency. 

WDL  is  the  prewarped  lower  frequency. 

! WDM  is  the  prewarped  center  frequency. 

i 

; WB  is  the  bandwidth,  WDU  - WDL. 

. T is  the  sampling  interval. 

The  Real  and  Imaginary  part  of  the  roots  of  the  filter  are  given  next. 

i I is  the  ith  stage.  I varies  from  1 to  N. 

I 

.•\q,  A^,  A2  are  the  Butterworth  band-stop  filter  numerator  coefficients, 
is  the  gain  factor. 

B^,  and  B^  are  the  Butterworth  band-stop  filter  denominator  coefficients. 
W is  the  frequency 

H is  the  overall  magnitude  of  the  digital  transfer  function 


HI 

is 

the 

magnitude 

of 

the 

digital 

transfer 

function 

(1st 

stage) 

H2 

is 

the 

magnitude 

of 

the 

digital 

transfer 

function 

(2nd 

stage ) 

H3 

is 

the 

magnitude 

of 

the 

digital 

transfer 

function 

(3rd 

stage) 

114 

is 

the 

magnitude 

of 

the 

digital 

transfer 

function 

(4th 

stage) 

See  Figure  4. 


b2 


This  is  an  example 
band-pass  filter  (0403) 
= 40  radians.  The  s 
increment  is  0.1  radian 


of  the  output  for  an  8th  order  Chebychev 
with  T = 0.002,  = 60  radians,  and 

tarting  frequency  is  0 radian,  the  frequency 
, and  the  ripple  is  0.1. 


WDU  = 60.07210  WDL  = 40.02135  WDM  = 49.02902  WB  = 20.05076 
T = 0.20000E-02 


A = 0.37642105 
A = 0.37642105 


B = 1.06850027 
B = 1.06850027 


K = 0.69553541 

Kg  = 0.28810020 


K = 0.28813942 
= 0.99524620 


THE  ROOTS  OF  THE  FILTER  ARE  GIVEN  BELOW 
REAL(l)  = -3.19528085  IMAGINARY(l)  = -44.97792290 

REAL(2)  = -3.77772492  IMAGINARY(2)  = -53.17662810 

REAL(3)  = -1.15829660  IMAGINARY(3)  = -40.10115290 

REAL(4)  = -1.73001696  IMAGINARY (4)  = -59.89456990 


THE  COEFFICIENTS  OF  EACH  DIGITAL  FILTER  SECOND  ORDER  SECTION  ARE 
GIVEN  BELOW 


FOR 

1=1 

= 

O.IOOOOOOOE+Ol 

= 

O.OOOOOOOOE+OO 

= 

-O.lOOOOOOOE+01 

= 

0.10395603E-01 

= 

-0.19792607E+01 

= 

0.98732565E+00 

FOR 

1 = 2 

= 

O.lOOOOOOOE+01 

s 

O.OOOOOOOOE+OO 

= 

-O.lOOOOOOOE+01 

= 

0.10375296E-01 

= 

-0.19737935E+01 

®2 

0.98504460E+00 

FOR 

1=3 

O.IOOOOOOOE+Ol 

O.OOOOOOOOE+OO 

-O.lOOOOOOOE+01 

= 

0.19406846E-01 

= 

-0.19889723E+01 

«2 

= 

0.99338493E+00 

FOR 

1 = 4 

A° 

O.lOOOOOOOE+01 

=: 

O.OOOOOOOOE+OO 

= 

-O.lOOOOOOOE+01 

0.19346636E-01 

= 

-0.19788675E+01 

«2 

= 

0.99312838E+00 

W 

0.0000  0. 
0.1000  0. 
0.2000  0. 
0.3000  0. 

0.4000  0. 

0.5000  0. 


II 

OOOOOE+00 

0. 

12493E-12 

0. 

19990E-11 

0. 

10121E-10 

0. 

31991E-10 

0. 

78116E-10 

0. 

HI 

OOOOOE+00 

51560E-03 

10312E-02 

15468E-02 

20625E-02 

25783E-02 


H2 

0. 00000 E+00 
0.36886E-03 
0.73773E-03 
0.11066E-02 
0.14755E-02 
0.18445E-02 


in 

O.OOOOOE+OO 
0.12106E-02 
0 24211E-02 
0.36318E-02 
0.48426E-02 
0.60536E-02 


114 

O.OOOOOE+OO 

0.54265E-03 

0.10853E-02 

0.16280E-02 

0.21707E-02 

0.27134E-02 
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WDU  is  the  prewarped  upper  frequency. 
WDL  is  the  prewarped  lower  frequency. 
WDM  is  the  prewarped  center  frequency. 
WB  is  the  bandwidth,  WDU  - WDL. 

T is  the  sampling  interval. 


K„  = 2Acos(0) 

O 

= A^cos^(6)  + B^sin“(9) 

The  Real  and  Imaginary  part  of  the  roots  of  the  filter  are  given  next. 

I is  the  ith  stage.  I varies  from  1 to  N. 

Aq,  A^,  A^  are  the  Chebychev  band-pass  filter  numerator  coefficients, 
is  the  gain  factor. 

and  B^  are  the  Chebychev  band-pass  filter  denominator  coefficients. 
W is  tile  frequency. 

H is  the  overall  magnitude  of  the  digital  transfer  function. 

HI  is  the  magnitude  of  the  digital  transfer  function  (1st  stage). 

H2  is  the  magnitude  of  the  digital  transfer  function  (2nd  stage). 

H3  is  the  magnitude  of  the  digital  transfer  function  (3rd  stage). 

H4  is  the  magnitude  of  the  digital  transfer  function  (4th  stage). 


See  Figure  5. 
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MAGNITUDE  VS  FREQUENCY 
FOR 

DIGITAE  TRANSFER  FUNCTION 


ORDER  BUTTERWORTH  BAND-PASS  FILTER 


N = 4 

Start  at  (1)  = 0 radian 
T = 0.002 

Steps  of  0.1  radian 
u)  = 60  radians 
(.1.  =40  radians 


ma(;nitui)K  vs  frequkncy 

FOR 

DIGITAI,  l'R.\NSFER  FUNCTION 
ORDER  BUTTERWORTH  BAND-STOP  FILTER 


N = 4 

Start  at  u)  = 0 radian 
T = 0.002 

Step.s  of  0.1  radian 
(.1  = BO  radian.s 

(1),  = 40  radians 


ma(;nitui)K  vs  i-rkouhncy 

I'OK 

DICITAI.  'I  R.\N.SKKK  Kl'NCTION 
OkDKR  CHRBYCHKV  HAND-I’ASS  KIl.TKR 


N = 4 

Start  at  . = 0 radian 

Ripple  = = 0.  100 

')'  = 0.002 

Steps  oi  0.1  radian 
iv  = hO  radians 
u).  = 40  radians 


MACNITUDK  VS  FRKQUl'.NCY 
FOR 

Dlon'Al.  TR.\NSFKR  FUNCTH)N 


ORDER  CHKBYCIIFV  RAND-STOF  FILTER 


iN  - ^ 

Start  at  m = 0 radian 
Ripple  = a = 0.100 
T = 0.002 

Steps  of  0.1  radian 
,.i  = 60  radians 

ai.  = 40  radians 
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INTRODUCTION 


This  report  contains  the  documentation  for  the  LPASS  program. 

It  consists  of  the  design  procedure  used,  a description  of  the 
program,  and  design  examples  using  the  program. 

The  purpose  of  the  LPASS  program  is  the  design  of  a maximally 
flat  Butterworth  or  an  equiripple  Chebychev  lowpass  digital  filter. 
Starting  with  an  analog  filter,  the  bilinear  Z transform  is  used  to 
find  an  equivalent  digital  filter.  The  user  enters  the  following 
parameters:  the  number  of  second  order  sections,  the  type  of  filter, 

the  sampling  interval,  the  -3db  cutoff  frequency,  the  starting  frequency 
and  the  frequency  increment.  If  a Chebychev  filter  is  being  designed, 
the  ripple  must  also  be  entered. 

The  program  calculates  the  digital  filter  coefficients  for  up  to 
three  second  order  sections  in  cascade.  The  program  is  designed  to 
calculate  up  to  a sixth  order  filter,  thus  the  filter  order  is  two 
times  the  number  of  cascaded  second  order  sections.  The  filter 
magnitude  response  is  generated  over  the  frequency  interval  specified 
by  the  input. 

The  LPASS  program,  written  in  Fortran  IV,  is  supplied  as  a card 
deck  with  this  report.  The  program  is  in  tlie  form  of  a subroutine 
and  can  be  used  as  is  by  a call  statement  from  the  m.iin  program. 

Data  may  be  input  via  cards  with  output  available  through  a line 
printer.  The  input/output  devices  may  be  altered  as  explained  in 
this  report.  Graphics  routines  may  easily  be  appended  to  t lu'  program. 


I. 


Design  Procedure 


A.  Preliminary  Discussion 

The  transfer  function  of  a second  order  digital  filter  in  the 
Z domain  is  given  by 


H(Z)  = 


K^(AqZ  +A^Z+A2) 
Z^+Bj^Z+R^ 


(1) 


where  the  A's  and  D's  are  the  coefficients  of  tlie  numerator  and 
denominator  respectively.  One  common  method  of  designing  a digital 
filter  is  to  start  with  an  analog  transfer  function  H(S)  and 
transform  it  to  the  digital  transfer  function  H(Z).  This 
program  will  calculate  the  scale  factor  and  the  coefficients 

Aq,  Aj^,  k^,  and  B,.  The  transformation  used  is  the  extendeu 

bilinear  Z transform  defined  as 


S ->■  — 

T ^z+r  ' 


(2) 


where  T is  the  sampling  interval.  Wtien  this  transform  is  employed  , 
the  desired  frequencies  must  first  be  prewarped  to  make  them  com- 
patible with  the  digital  filter.  The  prewarped  cutoff  frequency  is 
given  by 

? 10  T 

WDC  = - Tan  (-y-)  . (3) 

This  prewarping  is  done  by  the  program. 

B.  Butterworth  Low-Pass  Filter 

We  start  with  a normalized  second  order  low-pass  filter  in 
the  S plane. 


H(S) 


1 

2Scos0  + 1 


(A) 
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where  the  ai.gLe  0 is  in  degrees  (in  the  program).  0 may  be  found 
from  tiie  Butterworth  circle  and  the  relationship 


S = e 


ijTi(2m-l)/2n 


(5) 


where  n is  the  order  of  the  filter  and  m = 1,  2,  3,  ...  , n. 

Tills  relationship  is  determined  by  the  following  procedure.  By 
definition,  a filter  is  n*"^  order  Butterworth  low-pass  if  its  gain 
characteristic  is 

. 2 


(6) 


1 ^ (f) 

c 


2n 


where  a is  the  gain,  to^  is  the  desired  cutoff  frequency  and  n is 
the  order  of  the  filter.  Note  that  goes  to  zero  as  u goes 

to  infinity,  indicating  the  filter  does  attenuate  the  higher  frequencies. 
To  determine  its  efficiency  as  a low-pass  filter  we  calculate 

2n-l 


1 to 

d I an  c 

'l+(-^) 
to 

c 


2n' 


3/2 


(7) 


Thus 


_d 

dll) 


|H^(jo))  I 


= 0 


D=n 


(8) 


for  all  n and  hence  t'.ie  gain  characteristic  stays  flat  for  lo  close 
to  0.  Also 


2o)  /2 


(9) 


and  hence,  the  decline  rate  or  "roll-off"  of  the  gain  characteristic 
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i 


at  0)  = becomes  sharper  as  n increases.  In  other  words,  the 
approximation  to  the  ideal  low-pass  filter  Improves  for  larger  n. 

The  order  n is  chosen  according  to  desired  specifications.  The 
references  have  equations,  curves,  and  tables  that  select  n,  given 
the  specifications.  For  example,  page  227  of  Rabiner  and  Gold  gives 
an  equation  for  calculating  n when  the  transition  band  is  specified. 

In  the  design,  the  poles  for  the  full  frequency  response,  H(S), 
of  the  n*"^  order  Butterworth  filter  must  be  determined.  Tlie  pro- 
cedure is  as  follows: 


[H(S)H(-S)]  = 

b-JO) 


1 + 


2n 


S 1 + (^) 
u w = — iu 

c j c 


1 + 


2-in  - 


for  n even 


2 

w 

>-  c J 


1 + - 


for  n odd 


1 - 


2 

L c J 


Setting  the  denominators  equal  to  zero. 


5 

(10)  • 


U) 

c 


(±1) 


l/2n 


(11) 


Tlius,  the  pole  locations  are  the  2n  roots  of  ±1,  depending  on  whether 
the  order  is  odd  or  even.  These  roots  are  located  on  a circle  with 
radius  centered  at  the  origin  of  the  S plane  and  have  symmetry 
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with  respect  to  both  real  and  imaginary  axes.  For  n odd,  a pair  of 
roots  arc  on  the  real  axis  and  the  rest  are  separated  by  r/n  radians. 
For  n even,  a pair  of  roots  are  located  Ti/2n  radians  from  the  real 
axis  and  tlie  rest  are  again  separated  by  n/n  radians.  No  roots  are 
on  the  imaginary  axis  for  either  even  or  odd  n. 

Let  Pj,  . . . , p^^  be  the  roots.  From  the  symmetry  of  the  pole 
locations,  if  p^,  . . . , p^  are  the  roots  lying  in  the  right-half 


plane,  the  left-half  plane  roots  are 


• • • 


-p  . llic 
n 


magnitude-squared  function  can  then  be  written  as 
H„(S)11  (-S)  = 


a (-1)  01  2 
c 


(S+p^)  . . . (S-Pj)  . . . (S-Pj^) 


(12) 


To  be  stable,  must  have  all  its  poles  in  the  left-half  plane,  thus 


aoj 

”n^^^  = . . (S+p^) 


(13) 


The  program  is  written  with  unity  gain  at  DC,(oi=0),  therefore  a = 1. 

In  order  to  locate  the  poles  as  specified  above,  consider  the 
following  set  of  equations. 


±jn(2m-l) 

, = -e 

m = 1 , 2 , . . . 

, n;  for  n 

even 

-1  = 

k = 0,  1,  . . . 

, n;  for  n 

odd 

(14) 

.Substituting  equations  (14) 
1-S,  = _^ijTT(2m- 

(jO  * m 

into  equation  (11) 
, m = 1 , 2, 

yields 

• • • , n; 

for  n 

even 

c 

'oi  'ik 

, k = 0,  1,  , 

• • . . n ; 

for  n 

(15) 

odd 

c 

Equations  (15)  will  give  the 

pole  locations  as 

described 

abovi.  . 

Consider  the  form  of  equations  (15) 


S = -u)  e 
c 


o)^|-cos0  1 jsin9| 


(Ib) 
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From  this  relationship,  it  can  be  seen  that  the  magnitude  for  each  pole 

is  u)^,  regardless  of  the  angle,  and  thus  all  the  poles  lie  on  a circle 

with  radius  co  . 

c 

j As  an  example  consider  a second  order  filter,  n = 2. 

J 

f S,  ^ _^±jTT(2m-l)/4  m = 1,  2 

^co  ±m 
c 

) 

I “ /±135° 


e = 45° 


The  relationship  of  these  roots  about  the  circle  of  radius  w 

c 

is  illustrated  in  Figure  1.  The  angle  6 is  always  measured  from 
the  negative  real  axis. 

In  the  program,  only  the  angle  (s)  less  than  90°  are  considered 
so  that  poles  lie  in  the  left-half  plane  since  poles  in  the  left-half 
plane  are  stable.  Putting  6 = 45°  into  equation  (4)  yields  poles 
at  -0.707  ± j0.707.  These  locations  are  in  the  left-half  plane. 

In  the  program,  only  even  order  filters  are  considered. 

Below  are  the  values  of  0 for  1,  2,  and  3 second  order  sections 
in  cascade. 


Cascaded 
Sect  ions 

N 


1 

O 


Filter 

Order 

n 


Angle 

G 


4 

6 


45° 

22.5°,  67.5° 
75°,  45°,  15° 


Tliesc  calculated  angles  are  incorporated  in  llic  program  in  the  order 
, iti  ihove. 

7H 


F 


For  N second  order  sections  there  are  N 0's.  Only  one  specific 
0 is  used  per  stage,  because  each  stage  has  only  one  set  of  pole 
locations. 

The  following  is  the  procedure  to  derive  the  magnitude  of  the  i*”^ 
stage,  wliere  i varies  from  1 to  N. 

Given  the  normalized  second  order  low-pass  transfer  function 
equation  (4),  we  employ  the  low-pass  to  low-pass  transformation  for 
an  arbitrary  cutoff  frequency  given  by 

S - — (17) 

w 

c 

For  the  i*''^  stage,  equation  (4)  becomes 

2 

0) 

H.(S)  = ^ 2 (18) 

S +2So)  COS0.+CJ 

C 1C 

The  extended  bilinear  Z transform,  equation  (2),  is  used  to  get  to 
the  digital  domain.  Employing  equation  (2)  on  equation  (18)  and  sub- 
stituting WDC  for  yields 


"l<^>  ■ 43 


WDC^(Z^+2Z+1) 


(19) 


-y(Z"-2Z+l)4^(Z^-l)WDC  cos0.+WDC^(Z^+2Z+l) 

T^  1 1 

Putting  the  denominator  of  equation  (19)  in  monic  form  yields  the 
transfer  function  for  the  i*"^  stage  of  the  filter 

K (A  Z^+A  Z+A  ) 

H (Z)  = --y- . (20) 

z -pb^.z-fb^, 


Equation  (20)  is  the  same  as  equation  (1)  with  the  exception  of  the 
subscripts.  In  equation  (20) 
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1 


Aq  A2 


^ = 2 


= -|  + WDC  cos6^  + 


(21) 


WDC 

Hi  " G. 

1 


«li  = 


2WDC^  - yI 


«2i  = 


G. 

1 

-A  - % WDCcose.+WDC^ 
^2  T 1 


ST 

Letting  Z = e and  S = ju  and  taking  the  magnitude  of  H^(jco)  we 
have 


^(A  cos (2wT)+A. cos (wT)+A„)^+(AQSin (2uT)+A^sin (wT) ) 

|h_.  (jw)  I = 


J~(cos  (2wT)+Bj^^cos  (toT)+B2^)  ^+(Sin  (2wT)+B^^sin  (wT)  ) 


(22) 


This  magnitude  function  is  the  same  for  both  the  Butterworth  and  the 
Ghebychev  filters  wliere  i varies  from  1 to  N. 

C.  Ghebychev  Low-Pass  Filter 

The  advantage  of  the  Ghebychev  low-pass  filter  over  the 
Butterworth  low-pass  filter  is  tliat  the  transition  band  of  the 
response  at  frequencies  greater  ttian  is  sharper  for  the  Ghebychev 

low-pass  filter.  This  is  achieved  by  specifying  a small  percentage 
of  ripple  in  the  low-pass  region.  Tlie  amplitude  of  the  ripple  is 
specified  by  the  quantity  S (labeled  RIP  in  the  program).  Figures  (■> , 

7 and  8 illustrate  the  rippling  for  second,  fourth,  and  sixth  order  fil- 
ters, respectively.  The  poles  of  the  filter  are  found  on  an  ellipse 
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described  by  two  Butterworth  circles  of  radii  A and  B with  A<B. 
The  location  of  the  poles  on  the  ellipse  is  a function  of  the  ripple, 
6,  and  is  given  by  the  following  equation: 


B,  A 


4(( 


r;T.  -ia/2N 
+l+e  ) 


"+l+e 


(23) 


where 


e 


1 

(1-6F 


(24) 


B is  given  for  the  plus  sign  and  A for  the  minus  sign.  The  Chebychev 
ellipse  then  has  major  axis  B and  minor  axis  A.  The  location  of  the 
S plane  poles  on  the  ellipse  is  given  by 


Real  Part  = A cos9 
Imaginary  Part  = B sinO 


(25) 


The  0's  are  the  same  as  given  for  the  corresponding  order  Butterworth 
filter.  An  example  of  Chebychev  pole  locations  is  illustrated  in 
Figure  2.  For  A = 1/2  and  B = 1 in  a fourth  order  filter, 

0 = 22.5°  and  67.5°.  The  ;hebychev  pole  locations  are  determined 
from  equations  (23),  (24),  and  (25). 

The  analog  second  order  Chebychev  low-pass  filter  is 
K„a 

H(S)  = (26) 

S +KgS+K., 


where 


a 


(l+c^)^*^^ 


1/N 


(27) 


E is  calculated  from  equation  (24)  and  N is  tlie  number  of  second 
order  sections.  and  are  calculated  by 
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2ACOS0 


(28) 


''S  = 


2 2 2 2 

K2  = A cos  6 + B sin  0 . 


(29) 


The 


substitution  of  the  low-pass  to  low-pass  transformation  for  some 


jtoff  frequency  u , equation  (17),  into  equation  (26)  yields 

2 


H(S)  = -j 


K„aca 
2 c 


(30) 


S + SKqOj  + oj  K« 
8 c c 2 


Using  the  extended  bilinear  Z transform,  equation  (2),  and  substituting 


WDC  for  w we  have  for  any  section 
c 


H(Z) 


K^aWDC^CZ^+ZZ+l) 


;-(Z^-2Z+l)  + |KgWDC(Z  -1)  + WDC^K2(z'^+2Z+1) 


(31) 


th 

Collecting  terms  yields  the  following  for  the  i section 


Hi<Z)  . 


Z +B, ,Z+B„. 
li  2i 


(32) 


where 


Aq  ^2  ~ ^ 


A^  = 2 


Gi  = -|  + I WDC- Kg  + WDC^K^ 


aK2WDC 


'li 


(33) 


'11 


2WDC^K  - 

T 

G. 
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-y  - I WDC-Kg  + WDC^K2 


1 varies  from  1 to  N.  These  coefficients  are  used  to  find  |H^(jto)| 
given  in  (22). 

For  applications  where  a sharper  roll-off  is  required  the  Chebychev 
filters  are  used.  The  roll-off  increases  with  n for  any  fixed  e. 

For  fixed  n,  the  roll-off  decreases  as  e decreases.  For  small 
£ the  ripple  width,  6,  is  small,  see  equation  (23),  but  so  is  the 
roll-off.  For  larger  e the  roll-off  improves  but  the  ripple  width 
increases.  In  the  first  case  the  filter  will  be  good  at  DC  and 
low  frequencies,  unsatisfactory  at  high  frequencies.  The  converse  is 
true  in  the  second  case. 

The  above  observations  suggest  the  procedure  to  be  used  in 
selecting  a Chebychev  filter  to  match  a set  of  specifications.  The 
permissible  ripple  width  specifies  e.  With  e fixed,  select  n 
to  attain  the  required  roll-off. 


1 1 . Using  the  Program 

The  first  data  card  read  into  the  program  contains  the  number  of 

second  order  sections  to  be  cascaded,  N,  and  the  type  of  filter 

desired,  KN.  N is  equal  to  1,  2,  or  3,  which  corresponds  to 

the  2'^'^,  or  order  filter  respectively.  KN  = 1 yields  a 

Butterworth  filter,  wliile  KN  = 2 yields  a Chebychev  filter.  The 

format  on  the  N,  KN  card  is  212.  The  second  data  card  read  in 

is  the  sampling  interval  T in  F10.6  format.  Wfien  choosing  T,  1/T 

should  be  approximately  equal  to  ten  times  the  cutoff  frequency,  . 

The  tliird  data  card  contains  the  value  of  w in  F10.4  format. 

c 

For  the  Butterworth  low-pass  filter,  is  the  -3db  cutoff  frequency. 

9 1/2 

For  the  Chebychev  filter  the  magnitude  of  the  response  is  1/(1+l“) 

= 1-6  at  u)  = . 0)  is  in  radians.  0 is  the  ripple  factor. 

If  the  desired  filter  is  Chebychev,  i.e.,  KN  = 2,  the  next 

data  card  is  the  ripple  factor  (RIP)  in  F5.3  format.  The  filter 

response  for  all  even  order  Chebychev  low-pass  filters  passes  through 

2 1/2 

1/ (1+E  ) =1-6  for  w = 0 and  w . For  odd  order  filters,  the 

c 

2 1/2 

magnitude  is  1 for  u = 0 and  l/(l+e  ) =1-6  for  u = oi  . 

c 

This  program  produces  only  even  order  filters.  If  the  desired 
filter  is  Butterworth,  i.e.,  KN  = 1,  this  data  card  is  omitted  from 
the  data  deck. 

The  final  data  card  is  the  starting  frequency  (FREQl)  and  the 
frequency  increments  (BELT)  in  radians.  The  format  of  the  FREQl, 

BELT  card  is  2F10.4.  Betermlne  BELT  by  the  following: 


BELT  = 


final  frequency  - starting  frequence 


This  is  necessary  because  there  are  1024  frequency  data  points  calculated 


in  the  program.  Choose  FREQl  and  BELT  to  insure  that  calculated  values 
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will  Include  the  data  of  interest.  For  maximum  effiriency  of  the 


program,  DEI.T  should  be  a multiple  of  2 so  no  decimal  to  binary 
conversion  errors  are  Incurred. 

The  digital  filter  coefficients  are  computed  and  printed  out  or 
each  second  order  section.  The  full  filter  magnitude  response,  as  well 
as  each  section  magnitude  response,  is  printed  for  each  of  the  frequency 
Increments  specified.  When  N = 1,  the  section  magnitude  response  is 
the  full  filter  magnitude  response  and  is  only  printed  once. 

The  program  may  be  easily  modified  to  incorporate  a graphics 
display  of  the  magnitude  response.  There  is  a comment  card  in  the 
LPASS  program  indicating  where  the  graphics  subroutine  call  card  should 
be  inserted. 

The  program  is  written  with  input  obtained  via  device  4 and  output 
written  to  device  6.  These  numbers  should  be  assigned  to  the  appro- 
priate devices  prior  to  running  the  program. 

The  program  was  developed  on  the  PDP-11/20  with  a DOS/BATCH 
operating  system.  Trial  runs  frequently  used  a TTY  terminal  as 
well  as  a card  reader  for  input  (device  4) ; and  a TTY  terminal  as 
well  as  a line  printer  for  output  (device  6) . 

Double  precision  arithmetic  is  employed.  To  decrease  required 
memory  storage,  only  the  frequency  interval  values  and  the  full 
magnitude  response  are  saved.  The  section  magnitude  responses  are 
printed  out,  but  are  not  stored.  The  program  will  produce  approxi- 
mately 21  pages  of  output. 

Sliown  below  are  sample  deck  set-ups  for  the  Chebychev  and 
Butterworth  low-pass  filters. 
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Format  Example 


1 

212 

0302  (3  sections  Chebychev  low-pass) 

2 

F10.6 

0.001  (T  = 0.001) 

3 

F10.4 

100  (oj  = 100  radians) 
c 

4 

F5.3 

0.10  (Ripple  amplitude  = 0.10) 

5 

2F10.4 

70  0.06  (Start  at  w = 70.  Steps  of 
0.06  radians.  Will  finish  just 
past  w = 131  radians.) 

1 

212 

0201  (2  sections,  4*"^  order,  Butterworth) 

2 

F10.6 

0.005  (T  = 0.005) 

3 

F10.4 

20  (w  = 20  radians) 
c 

4 

2F10.4 

0 0.04  (Start  at  co  = 0,  finish  just 
past  (jj  = 40  radians  in  steps  of 

0.04  radians.) 


The  following  pages  contain  annotated  examples  of  output  data. 
This  is  an  example  of  the  output  for  a 4*"''  ofder  Butterworth 
low-pass  filter  with  T = 0.005  and  = 20  radians.  The  starting 
frequency  is  0 radians  and  the  frequency  increment  is  0.04  radian. 


WDC  = 20.01668  WC  = 20.00000  T = 0.50000E-02 


FOR  I = 1 Aq  = 


FOR  1=2  Aq  = 


O.lOOOOOOOE+01  A^ 
O.lOOOOOOOE+01 
0.18219614E+01  B2 

O.lOOOOOOOE+01  A^ 
O.lOOOOOOOE+01 
0.19167786E+01  B2 


= 0.20000000E+01 
= 0.22869799E-02 
= 0.83110937E+00 

= 0.20000000E+01 
= 0.24059972C-02 
= 0.92640257E+00 


W 


H 


(12 


0.0000 

0.0400 

0.0800 

0.1200 

0.1600 

0.2000 

0.2400 


O.lOOOOE+01 
O.lOOOOE+01 
0.99999E+00 
0.1 0000 E+01 
O.lOOOOE+01 
0. lOOOOE+01 
O.lOOOOE+01 


O.lOOOOE+01 
O.lOOOOE+01 
0.99999E+00 
0. 99997 F.+OO 
0. 99995E+00 
0.99993E+00 
0.99990E+00 


O.lOOOOE+01 
O.lOOOOE+01 
0. lOOOOE+01 
0.  lOOOOF.+Ol 
O.lOOOOE+01 
O.lOOOOE+01 
0. lOOOlE+01 


1 


I is  the  stage.  I varies  from  1 to  N. 

WDC  is  the  prewarped  cutoff  frequency. 

WC  is  the  cutoff  frequency. 

T is  the  sampling  interval. 

Aq,  A^,  and  .\2  are  the  low-pass  filter  numerator  coefficients, 
and  B2  are  the  low-pass  filter  denominator  coefficients, 
is  the  gain  factor. 

W is  the  frequency. 

H is  the  overall  magnitude  of  the  digital  transfer  function. 

s t 

HI  is  the  magnitude  of  the  digital  transfer  function  (1  stage). 

H2  is  the  magnitude  of  the  digital  transfer  function  (2’^'^  stage)  . 

H = H1*H2. 

See  Figure  4. 

This  is  an  example  of  the  output  for  a order  Chebychev 
low-pass  filter  (three  second  order  stages  cascaded)  with  T = 0.005 
and  = 20  radians.  The  starting  frequency  is  0 and  the  frequency 
Increment  is  0.04  radian.  The  ripple  is  equal  to  0.100. 

WDC  = 20.01668  WC  = 20.00000  T = 0.50000E-02 

A = 0.24783947  B = 1.03025433  K„  = 0.12829114  K = 0.99443709 
A = 0.24783947  B = 1.03025453  k“  = 0.35049793  = 0.56142438 

A = 0.24783947  B = 1.03025453  = 0.47878908  = 0.12841170 

o L 

FOK  I = 1 Ap  = O.lOOOOOOOE+01  A^^  = 0. 20000000E+01 

= 0. lOOOOOOOE+01  = 0.23830688E-02 

B^  = -0.19774006E+01  B2  = 0. 98727357E+00 

WDC2  = 0.40066761E+03  G(I)  = 0. 16142562E+06  A = 0. 96548939E+00 

FOR  1=2  Aq  = 0. lOOOOOOOE+01  A^  = 0. 20000000E+01 

A2  = O.lOOOOOOOE+01  Kj  = 0.13321469E-02 

Bj^  = -0.  1960054  lE+Ol  B2  = 0. 96557320E+00 
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WDC2 


0.40066761E+03 


G(I)  = 0.16303127E+06 


A = 0.96548939E+00 


FOR  I 


Aq  = O.lOOOOOOOE+01 
A^  = O.lOOOOOOOE+01 
= -0.19519613E+01 


A^  = 0.20000000E+01 
= 0. 30310788E-03 
= 0.95321709E+00 


WDo2  = 0.40066761E+03  G(I)  = 0. 16388496E+06  A = 0. 96548939E+00 


w 

H 

Hi 

H2 

H3 

0.0000 

0.89998E+00 

0.96549E+00 

0.96549E+00 

0.96547E+00 

0.0400 

0.89999E+00 

0.96549E+00 

0.96549E+00 

0.96548E+00 

0.0800 

0.90001E+00 

0.96550E+00 

0.96551E+00 

0.96547E+00 

0.1200 

0.90009E+00 

0.  96552 E+00 

0.96554E+00 

0.96550E+00 

0.1600 

0.90016E+00 

0.96555E+00 

0.96558E+00 

0.96551E+00 

0.2000 

0.90028E+00 

0.96558E+00 

0.96564E+00 

0.96555E+00 

0.2400 

0.90042E+00 

0.96563E+00 

0.96571E+00 

0.96559E+00 

WDC  is  the  prewarped  cutoff  frequency. 

WC  is  the  cutoff  frequency. 

T is  the  sampling  interval. 

D A 1 ///  -2.,^  -1,1/2N^, -1.-1/2N, 

B,  A = — ((^e  +1+C  ) ±(\le  +l+e  ) ) 

I is  the  i*"^  stage,  I varies  from  1 to  N. 

Kg  = 2Acos6. 

2 2 2 
K2  = A cos  0+B  sin  0. 

Aq,  A^,  and  A^  are  the  low-pass  filter  numerator  coefficients. 
Bj^  and  B2  are  the  low-pass  filter  denominator  coefficients. 

is  the  gain  factor. 

WDC2  = (WDC)^. 

C(I)  = -^  + + (WDC)^K2. 


The  A following  G(I)  is  a 


Li  + c"J 


1/2N 


W is  the  frequency. 

H is  the  overall  magnitude  of  the  digital  transfer  function. 
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s t 

HI  is  the  magnitude  of  the  digital  transfer  function  (1  stage). 

H2  is  the  magnitude  of  the  digital  transfer  function  (2^^  stage). 

rd 

H3  is  the  magnitude  of  the  digital  transfer  function  (3  stage). 
H = H1*H2*H3. 


See  Figure  8. 
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ABSTRACT 

Given  an  Nth  order  discrete  system  the  autocovariance  sequence  is 
formulated  as  a set  of  recursive  and  nonrecurslve  difference  equations. 
Application  of  the  final  value  theorm  for  a z transform  then  permits  the 
first  N terms  of  the  steady  state  autocovariance  sequence  to  be  computed 
via  one  matrix  inversion.  The  remaining  terms  are  calculated  using  a 
simple  recursive  relationship. 


99 


INTRODUCTION 


In  the  design  of  discrete  systems  it  is  often  necessary  to  compute 
the  variance  or  the  autocovariance  due  to  input  and/or  quantization  noise. 

In  many  instances  these  noise  sources  can  be  represented  as  independent 
zero  mean  sequences  [1]  which  allows  the  steady  state  auto- 
covariance  sequence  {R[kT]}  to  be  represented  by  the  discrete  Wiener- 
Khinchine  relationship  [1,2] 

R[kT]  = ~ 0 H[z]  H[z"^]  z*^"^  dz  . (1) 

2 

In  (1)  a is  the  variance  of  the  noise  and  H[z]  is  the  discrete 
transfer  function  from  the  noise  source  to  the  system  output.  For  a noise 
sequence  starting  at  t=0,  the  transient  autocovariance  sequence 
{R  [kT]}  time  t=nT  is  given  by 

n-k 

R IkT]  =0  y h[(l+k)T]  hllT]  (2) 

i=0 

where  h[nT]  is  the  inverse  z transform  of  l'.[z].  As  m becomes  large, 

each  terra  R [kT]  in  the  autocovariance  sequence  converges  to  R[kT] 
m 

given  by  (1) . 

From  a computational  point  of  view,  (1)  can  be  evaluated  using  the 
residue  theorm  from  complex  variables.  If  only  the  steady  state  variance 
is  of  interest.  Jury  [2]  has  derived  an  expression  for  R[OT]  that  is 
the  ratio  of  two  determlnents  whose  terms  are  simple  functions  of  the 
transfer  function  coefficients.  The  use  of  (2)  requires  the  evaluation 
of  the  weighting  sequence  and  a series  summation  for  all  desired  values 
of  k . 
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In  this  paper  the  autocovariance  sequence  is  formulated  as  a set 


of  difference  equations.  If  only  the  steady  state  autocovariance  is  of 
interest,  the  use  of  z transforms  leads  to  a set  of  simultaneous  equations 
that  can  be  easily  evaluated  via  one  matrix  inversion.  Remaining  terms 
in  the  autovariance  sequence  are  found  using  a simple  nonrecursive 
relationship . 

DIFFERENCE  EQUATIONS  FOR  THE  AUTOCOVARIANCE 
First  consider  an  Nth  order  discrete  transfer  if^nctlon  of  the  form 

W 

H[z]  = . (3) 

1 + I b z'^ 
j=l  ^ 

For  an  Independent  zero  mean  input  noise  sequence  {w[nT]}  with  variance 
2 

ow  that  starts  at  t=0,  the  difference  equation  describing  the  system  is 

N 

y[nT]  = - I b.  y[(n-j)T]  + w[nT]  . (4) 

j=l 

The  mean  value  of  y[nT]  is  zero  for  all  n and  terms  in  the  autocovariance 
sequence  at  time  t=nT  are  now  defined  as 


R^[kT]  = E{y[nT]  y[(n-k)T]}  k=0,l,2,...n  . (5) 

From  (4)  and  (5)  the  variance  of  y[nT]  is  given  by  the  quadratic  form 


R [OT]  = b*^  Q b + ow^ 
n n 


(6) 


where  is  a time  varying  symmetric  autocovariance  matrix  with  terms 


“^ij  ° *^n-min(l,J) 

= 0 
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n-min(i,j)  ^ 0 
n-mln(l,j)  < 0 


(7) 


The  vector 


The  vector  is  defined  as  the  row  matrix  [-b  -b  ...-b  ] . 

12  n 

Since  w[nT]  is  taken  from  an  independent  noise  sequence,  y[(n-k)T] 
and  w[nT]  are  uncorrelated  for  k^l  and  the  next  N terms  in  the  auto- 
covariance sequence  are  found  by  multiplying  both  sides  of  (4)  by 
y[(n-k)T],  k=l,2,...N  and  taking  the  expectation.  This  yields  a set  of 
equations  described  in  vector  form  as 

R = Q b (8) 

n 'n  ' 

where  R ^ is  the  row  matrix  |R  [T]  R [2T]  ...  R [NTll  . 

n U n n ■'  n J 

The  N+1  equations  described  by  (6)  and  (8)  consist  of  a set  of  N/2 
coupled  recursive  difference  equations  for  N even  and  (N+l)/2  coupled 
recursive  equations  for  N odd.  The  remaining  equations  are  nonrecursive 
difference  equations.  For  k>N  the  R^[kT]  terms  are  given  by 

N 

R [kT]  = - [ b.  R^  . [(k-l)T]  . (9) 

n i n— 1 

A closed  form  steady  state  solution  is  derived  by  first  defining  the 

z transform  of  the  autocovariance  term  R .[kT]  as 

n-i" 

Z {R  , [kT] } = z~^  R [kT]  (10) 

n-i  z 

where  the  subscript  (n-i)  represents  the  discrete  time  (n-l)T. 

Taking  the  z transform  of  (6)  and  (8)  now  yields  the  expression  for  the 
z transform  of  the  variance 


R^[0T]  = b*^  Q[z]  b + ow^ 


and  the  next  N terms  of  the  autocovariance  sequence 


(11) 


(12) 
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In  (11)  the  matrix  Q[z]  has  terms 


[zi  - R^nniu 

where  1 1- j ! takes  on  values  0 , 1 , 2 , . . ,N-1 . Expanding  the  quadratic  form 
given  by  (11)  and  applying  the  final  value  theorm  to  this  and  the  first 
N-1  equations  given  by  (12)  now  yields  a set  of  N simultaneous  equations 
whose  solution  forms  the  first  N terms  of  the  steady  state  autovariance 
sequence.  In  matrix  form  we  have 


RtOT] 


R[(N-1)T]| 


where  M is  the  matrix 


1 - I I>1  -2  I !>  >.  -2  I b b 

i=l  1=1  1=1  ^ ^ 


M=  (b^) 


(l+b^) 
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The  remaining  terms  in  the  autocovariance  sequence  are  given  by 


R[kT]  = - )'  b,  R[(k-i)T]  k > N. 

^ j 


i=l 


(16) 


For  a system  transfer  function  of  the  form 
V-1 


H[z]  = 


I a 2 

i=0 


1 + I N 


_ = 

N -1 

-j  D[z 


(17) 


j=0 


J 


the  transfer  function  is  separated  by  defining  an  intermediate  operation 


r -1,  W[z"^] 

x[z  ] = - 

D[z 


(18) 


so  that  the  response  is  given  by 


y [z  = Nf  z x[z 


(19) 


The  autovariance  sequence  for  x^  is  found  using  the  proceeding  results. 

The  terms  in  the  output  autocovariance  sequence  are  given  by 

E{y[nT]  y[(n-k)T]}=  a (20) 

where  a*"  is  the  row  matrix  [a.  a^  . . . a„  , ] and  C is  the  appropriate 

0 1 V-1  n rr 

V by  V time  varying  matrix.  In  steady  state  is  a constant  matrix 

with  terms 


C^.  = R[(k-j+i)T]  = Lim  E{x[nT]  x [ (n- (k-j+i) )T] } . (21) 

1 n-«° 

When  computing  the  steady  state  autocovariance  for  K < V certain  terms 
in  the  C matrix  are  found  from  the  relationship  R[kT]  = R[-kT].  Thus, 


I 

J 

'i 
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r 


1 


I 

i 


a fixed  nonrecursive  expression  exists  only  for  K > V. 


EXAMPLE 


To  illustrate  the  procedure,  we  consider  a low  pass  fourth-order 
Butterworth  filter  with  a -3db  point  at  w^  = 10.  The  filter  was 
designed  using  the  bilinear  z transform  with  T=0.01  and  has  the  general 
transfer  function 


where 


and 


H[z] 


K[z  + 1] 


4 3 2 

z + b-z  + b-z  + b_z  + b, 
1 z J 4 


N[z] 

n[z] 


k = 4.69832343  x lO"^, 

b^  = -2.53346973, 

b^  = 2.65559567, 

b^  = -1.28757608, 

b,  = 0.24062331. 

4 


(22) 


For  a direct  form  of  implementation  [1]  with  nine  product  rounding 
errors  the  difference  equation  for  the  output  error  is 

4 

e[nT]  = p[nT]  - ^ b e[(n-j)T]  . (23) 

j=l  ^ 

2 

The  variance  of  p[nT]  is  9q  /12.  The  variance  of  e[nT]  was  computed 
as  a function  of  n using  (6)  and  (8)  and  is  plotted  in  Fig.  1.  This  was 
checked  using  the  form  given  by  (2)  for  k=0  and  n=“>.  The  steady  state 
autocovariance  for  the  filter  is  shown  in  Fig.  2.  Here  the  M ^ matrix 
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from  (14)  is  given  by 


64.818991 


57.510375 


38.381329 


-213.413104 


-188.201279 


-122.830831 


129.627082 


115.638784 


79.495504 


-27.676674 


-24.832401 


-17.239467 


. (24) 


14.134665 


- 40.901684 


33.388940 


- 6.391421 


The  multiplication  shown  in  (14)  gives  the  first  four  terms  of  the  sequence 
with  the  remaining  terms  given  by  (16). 

For  two  second  order  sections  in  series  given  by 

H [z]  = (5.78776100)10"^  (z-H)^  ^ ^1^^^ 

^ - 1.07350061  z 0.30805006  ^1^^^ 


r_i  (7.99359506)10"^(z-t-l)^  ^2^^^ 

H Izj  - — ..  (2f 

z - 1.45996913  z + 0.77971293  2*^-' 

the  implementation  of  each  section  in  direct  form  requires  the  following  two 


calculations  to  compute  the  autovariance  sequence  at  the  output  of  H2[z]. 

First  the  autocovariance  due  to  the  noise  in  section  one 

is  found  using  the  coefficients  in  the  transfer  function 


HijIzJ  . 


D^[z]  D2[z] 


Thus  the  autovariance  sequence  shown  in  Fig.  2 is  modified  by  the  zeros 
of  H2[z].  From  (20)  this  Is  given  by 


R^2tkT]  = a C a . 


The  autocovariance  R22lkT]  due  to  the  noise  in  section  two  using  the 
coefficients  from 


^22^^^  " D^[z] 


(29) 


The  steady  state  autovariance  of  the  filter  is  now  given  by 


Ry[kT]  = R^2f^T]  + R22[1^]  (30) 

since  the  product  rounding  errors  are  assumed  to  be  uncorrelated  [!]• 

A plot  of  R^[kT]  is  shown  in  Fig.  3 as  a function  of  k.  This  Illustrates 
the  well-known  result  that  implementation  of  high  order  filters  as  a 
series  or  parallel  combination  of  first  and  second  order  sections  can 
reduce  the  effect  of  product  rounding  errors. 


CONCLUSIONS 

A method  for  computing  the  autocovariance  at  the  output  of  a digital 
filter  has  been  presented.  For  the  steady  state  case,  the  inversion  of 
an  N by  N matrix  is  required  where  N is  the  order  of  the  system  from  the 
noise  source  to  the  filter  output.  The  remaining  terms  are  computed 
using  simple  nonrecurslve  expressions  involving  the  filter  coefficients. 

The  procedure  is  suitable  for  programming  on  a computer  or  on  one  of  the 
programmable  scientific  calculators  that  now  exist  in  many  design  facilities. 


V plotted  as  a function  of  n for 
n ^ 

eiTiented  In  direct  form. 


Steady  state  autocovariance  f 
filter  Implemented  as  two  sec 
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ABSTRACT 

The  variance  of  error  at  the  output  of  a digital  control  system 
due  to  input  data  quantization  and  product  rounding  in  a digital  con- 
troller is  derived.  Two  forms  of  controller  implementations  are 
considered.  Then  an  expression  for  the  word  length  that  yields  a 
specified  variance  of  error  is  determined.  Scaling  is  then  discussed 
and  an  example  consisting  of  a discrete  Integral  plus  proportional 
controller  is  presented. 
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INTRODUCTION 


The  use  of  discrete  data  processing  algorithms  as  an  Integral  part 
of  a control  system  is  well  known.  Books  by  Tou  [1],  Monroe  [2]  and 
Kuo  [3]  describe  many  of  the  design  concepts  for  these  systems  called 
sampled  data  or  digital  control  systems.  In  a somewhat  Independent 
effort  design  and  implementation  procedures  have  been  developed  for 
digital  filters.  A collection  of  papers  In  this  area  has  recently  been 
edited  by  Rabiner  and  Rader  [A].  A tutorial  paper  that  provides  an 
overview  of  errors  in  the  Implementation  of  digital  filters  has  been 
written  by  Liu  [5]. 

The  difference  between  a digital  filter  and  a digital  control  system 
is  in  the  system  structure.  A digital  filter  Is  inherently  an  open  loop 
signal  processing  component  that  is  described  entirely  by  a difference 
equation.  In  a digital  control  system,  the  discrete  control  algorithm, 
which  is  a digital  filter,  is  a component  in  a feedback  loop  that  also 
usually  contains  analog  elements.  Analysis  of  the  mixed  system  then 
requires  the  use  of  artificial  sampling  elements  within  the  system. 

Currently,  discrete  control  algorithms  are  often  implemented  using 
a digital  computer  that  is  programmed  to  handle  one  or  more  control  paths. 
However,  for  some  applications,  the  real-time  speed  requirements  can 
tax  the  capacity  of  a computer.  Furthermore,  the  real-time  programming 
of  many  parallel  control  algorithms  is  costly  and  difficult  to  manage. 

As  an  alternative,  the  recent  development  of  small  mlcroprograramable 
processors  now  permits  a different  design  philosophy  for  digital  control 
systems.  The  use  of  these  processors  to  implement  discrete  control  algorithms 


i 
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as  Independent  units  or  under  control  of  a central  computer  leads  to  a 
distributed  computing  structure  that  Is  well  suited  for  use  In  a variety 
of  control  systems.  In  a typical  system,  each  mlcroprogrammable  processor 
is  Initially  programmed  from  a central  computer  to  perform  a system  function. 
If  operating  conditions  change  the  processors  can  be  reprogrammed  to 
meet  new  requirements. 

In  the  Implementation  of  any  discrete  control  algorithm,  the  designer 
should  consider  the  effect  of  errors  due  to  finite  length  arithmetic.  The 
errors  that  are  Introduced  are  caused  by  the  quantization  of  the  Input 
data,  quantization  of  the  algorithm  coefficients  and  the  quantization  of 
products.  When  floating  point  hardware  Is  used,  sums  and  differences 
also  Introduce  errors.  For  satisfactory  system  operation,  the  algorithm 
should  be  Implemented  In  a form  that  minimizes  the  effect  of  these 
errors  and  the  word  length  needed  for  the  imnlementatlon  should  be  determined. 

Quantization  has  been  studied  by  a number  of  researchers.  Bennett  [6] 
and  Widrow  [7]  have  shown  that  for  a rounding  form  of  quantization,  the 
errors  can  be  modeled  as  independent  random  variables  with  mean  zero 
2 

and  variance  q /12.  Here  q represents  the  value  of  the  least 
significant  bit  In  the  computer  word.  This  model  has  been  experimentally 
verified  for  Input  signals  that  transverse  several  quantization  levels 
from  sample  to  sample.  The  truncation  form  of  quantization  Is  similar  to 
rounding  except  that  the  mean  value  for  each  error  Is  q/2. 

A method  for  finding  the  maximum  upper  bound  on  the  error  In  a 
digital  control  system  due  to  quantization  was  Introduced  by  Bertram  [8]. 
Slaughter  [9]  determined  an  upper  bound  on  the  steady-state  system  error 
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by  replacing  each  error  source  by  a step  function  equal  to  the  maximum 
value  of  the  error  source.  The  steady  state  error  is  then  evaluated  by 
I applying  the  final  value  theorem  from  z-transforras . Knowles  and  Edwards 

i 

i [10,  11]  assumed  random  quantization  errors  and  developed  statistical 

i 

1 upper  bounds  for  the  error  at  the  output  of  the  system. 

I 

, In  this  paper,  the  variance  of  the  error  at  the  system  output  due 

j to  rounding  of  the  input  data  and  the  Intermediate  products  in  a digital 

j control  algorithm  is  evaluated  for  two  common  algorithm  Implementations. 

A unity  feedback  is  assumed,  however,  the  procedure  can  be  easily 
extended.  Fixed  point  arithmetic  is  also  assumed  since  this  form  of 
arithmetic  is  common  in  small  microprocessors.  Because  scaling  is  also 
a consideration,  scaling  procedures  are  described  and  their  effect  is 
discussed.  Finally,  an  example  consisting  of  a discrete  form  of  an  Integral 
plus  proportional  controller  is  investigated  and  the  word  length  needed 
for  a given  variance  is  compared  to  that  obtained  using  methods  described 
in  reference  [9] . 

EVALUATION  OF  THE  ERROR  VARIANCE  IN 
DIGITAL  CONTROL  SYSTEMS 

A block  diagram  for  a unity  feedback  digital  control  system  is  shown 
in  Fig.  1.  Here  the  input  signal  x(t)  and  the  return  signal  from  the 
I plant  are  assumed  to  be  analog.  These  signals  are  converted  into  digital 

form  and  subtracted  to  form  the  discrete  error  signal  e[nT].  This 
signal  serves  as  the  input  to  the  digital  controller  with  transfer  function 
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Hill  IM  hm‘| 


C(z)  = 


1 + I b z 


Cp(z)  • 


Errors  In  the  Implementation  of  (1)  are  due  to  rounding  of  the  analog 

data  by  the  analog-to-d igital  converter  and  the  rounding  of  products. 

Coefficient  quantization  errors  are  deterministic  and  under  control  of 

the  designer.  These  errors  will  not  be  considered  here.  For  simpllcitv, 

the  least  significant  quantization  level  will  be  taken  as  + q/2 

for  all  errors.  This  is  not  a serious  restriction  and  different  quantization 

levels  can  be  used.  Each  rounding  error  is  also  assumed  to  be  an  Independent 

2 

random  variable  with  mean  zero  and  variance  q /12. 


Quantization  of  the  Input  Data 

For  the  system  shown  in  Fig.  1,  two  input  data  quantization  errors 
occur.  These  errors  are  represented  by  noise  sources  shown  in  the 
block  diagram  of  Fig.  2.  The  error  at  the  system  output  due  to  the  two 
errors  at  discrete  times  t=nT  is  given  bv 


e tnT]  = I h[kT]  (e  [(n-k)T]  - e [(n-k)T]} 

^ k=0  ^ 


where  h[kT]  is  the  system  weighting  sequence  given  bv  the  Inverse 
z transform  of  the  system  transfer  function 


‘ 1 + C[z]  G[z] 


In  (3),  G[z]  is  the  z transform  equivalent  of  the  data  reconstruction 


element  and  the  fixed  plant.  If  Cj^[kT]  and  62[kT]  are  independent 
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random  variables,  the  variance  at  the  system  output  at  time  t=nT  is 


Var{e  [nX]}=^  I h^[kT]  . (4) 

^ k=0 

The  expression  given  by  (4)  is  monotonically  increasing  as  n Increases  and 
approaches  steady  state  as  n approaches  infinity.  The  steady  state 
variance  is  also  given  by  the  discrete  Wlener-Khinchine  relationship 

Var  ^ H[z]  H[z  z ^ dz^  . (5) 

The  evaluation  of  the  steady  state  variance  and/or  the  steady-state  auto- 
variance  sequence  is  easily  accomplished  using  the  procedure  developed 
in  Appendix  1. 

Product  Rounding  Errors 

Product  rounding  errors  occur  within  the  digital  control  algorithm 
and  their  effect  on  the  system  depends  on  the  form  of  implementation. 

Typical  forms  are  the  direct  form  one  and  direct  form  two.  To  minimize 
product  rounding  errors,  high  order  algorithms  are  usually  implemented 
as  a cascade  or  parallel  connection  of  first  and  second  order  sections  [5]. 
The  error  analysis  for  these  cases  follows  the  general  results  presented 
in  this  paper  and  will  not  be  considered  further. 

The  block  diagram  for  a third  order  algorithm  Implemented  in  direct 
form  one  is  shown  in  Fig.  3 where  each  multiply  operation  consists  of  a 
rounded  product  plus  a noise  component.  For  this  form  of  implementation, 
the  transfer  function  from  each  error  source  to  the  controller  outout  is 


given  by 


For  an  Nth  order  algorithm  described  by  (1)  the  resulting  error  at 
the  output  of  the  control  system  for  a direct  form  one  implementation 


is  given  by 


where 


e ,[nT]  = [ h [kT]  p[(n-k)T] 

^ k=0  ^ 


p[(n-k)T]  = I e [(n-k)T]  + I 6 [(n-k)T] 

k=0  j = l J 


and  hj^[kX]  is  the  Inverse  z transform  of 

Cj^[z]  G[z] 

" 1 + C[z]  G[z]  • 


The  variance  of  the  error  for  Independent  product  rounding  errors  is 

0 

Var{e  _[nT]}  = (M+R)  ^ V h,  [kT]  . 

k=0  ^ 


The  variance  of  the  total  error  due  to  input  data  quantization  and 
product  rounding  for  the  digital  control  algorithm  Implemented  in  direct 
form  one  is  now  given  by 

2 r n N I 

Var{e  [bT]}=  ^ ^2  [ h [kT]+  (M+R)  [ h [kT]>  . (11) 

^ I k=0  k=0  ^ J 


The  other  common  form  of  Implementation  for  a digital  controller  is 
the  direct  form  two  shown  in  Fig.  4 for  a third  order  system.  Here,  the 
errors  generated  by  multiplying  with  the  denominator  coefficients  can  be 
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considered  as  input  errors  to  the  digital  controller.  The  errors  due  to 
the  a^  coefficients  are  simply  additive  at  the  controller  output.  The 
resulting  variance  of  error  at  the  system  output  due  to  input  data 
quantization  and  product  rounding  for  a direct  form  two  implementation  is 

2rn  '^2  ''2I 

Var{e  [nT]}  = ^2  [ h [kT]  + R I h^[kT]+  (M)  I h,[kT]i  (12) 

^ I k=0  k=0  k=0  ^ } 


where  h2[kT]  is  the  Inverse  z transform  of 


H2[z] 


C.[z] 

1 + C[z]  G[z] 


(13) 


WORD  LENGTH  REOUIREMENTS  USING  THE 
VARIANCE  OF  ERROR  AT  THE  SYSTEM  OUTPUT 

For  a given  analog  signal  range  of  + A,  the  level  of  least  quantiza- 
tion referred  to  the  algorithm  input  is  given  by  the  inequality 


In  (14),  N represents  the  number  of  magnitude  bits  that  are  combined  with 
a sign  bit  within  the  computer.  The  Inequality  is  used  because  the  least 
level  of  quantization  is  often  chosen  so  the  range  of  the  scaled  binarv 
variables  slightly  exceeds  the  value  of  A.  For  example,  with  a q of 
10  mv  and  a ten  bit  and  sign  word  length,  the  maximum  value  of  the 
binary  variables  in  terms  of  signal  units  is  10.23  volts.  Carefully  note 
that  q is  given  in  terms  of  the  analog  range  at  the  system  input. 
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Given  a specification  on  the  variance  of  error  at  the  system  output, 


the  word  length  that  allows  the  specif Icatioit  to  be  met  is  found  by 


substituting  (14)  into  (11)  or  (12)  and  solving  the  resulting  inequality 

2 

for  N.  For  a specified  steady  state  variance  , the  word  length 

requirement  to  meet  or  exceed  the  specification  is  found  by  substituting 

(14)  into  (11)  or  (12)  with  Var{e  „,[“]}  = a ^ or 

yDl  s 

2 

Var{e  = a . This  gives  the  word  length  for  a direct  form  one 

yuZ  s 

implementation  of  the  digital  control  as 
N'di  > log 


^ ^ I2  y h^[kX]  + (M+R)  I h ^[kT]|)‘'^+ll  . 
I\l2  a * k=0  k=0  ^ / J 


(15) 


For  the  direct  form  two  we  have 


I h^[kT]  + R I h^[kT]  + 
k=0  k=0 


(M) 


00 


k=o 


h2^[kT]|V'^+l 


(Ifi) 


In  (15)  and  (16)  and  ^re  the  smallest  Integers  that  satisfy 

the  inequalities.  In  addition  a sign  bit  is  normally  required. 


SCALING 

Along  with  word  length  considerations  for  input  data  quantization 
and  product  rounding,  the  designer  of  a digital  control  algorithm  must 
consider  scaling.  For  the  algorithm  to  operate  properly,  the  output  r(nT) 
as  shown  in  Fig.  1,  cannot  overflow  for  any  given  input  signal  x(t). 

If  this  happens,  the  control  loop  is  effectively  open  and  the  system 
will  not  function  in  a desired  manner. 
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For  the  common  case  of  fixed  point  two's  complement  arithmetic, 
Jackson,  Kaiser  and  McDonald  [11]  have  shown  that  intermediate  overflow 


at  the  output  of  any  multiplier  or  adder  will  not  alter  the  algorithm 
performance  as  long  as  the  algorithm  response  does  not  overflow.  This 
means  that  for  the  filter  form  shown  in  Fig.  3 the  response  r(nT)  is 
the  only  variable  that  must  be  scaled.  For  the  direct  form  two,  both 
r(nT)  and  r^(nX)  must  be  scaled.  Scaling  for  these  Implementations 
has  been  investigated  by  Jackson  [12]  who  used  Holders  inequality  to 
establish  steady  state  bounds  for  various  digital  filter  forcing  functions. 
However,  these  bounds  do  not  Include  the  effect  of  transients. 

In  practice,  analytical  methods  for  determining  bounds  on  r(nT) 
often  gives  pessimistic  results  and  the  use  of  simulation  for  a given 
class  of  input  signals  can  be  useful.  Once  the  maximum  value  or  a bound 
for  r(nT)  and/or  rj^(nT)  is  determined  the  actual  scaling  is  not 
i.  difficult.  For  the  direct  form  one  implementation  shown  in  Fig.  3,  each 

a^  coefficient  is  reduced  in  magnitude  by  a common  factor.  This 
effectively  reduces  the  gain  of  the  algorithm.  For  a direct  form  two 
implementation,  the  input  signal  e(nT)  is  multiplied  by  a scale  factor 
so  that  the  values  of  r^(nT)  do  not  overflow.  The  a^  coefficients  are 
reduced  in  magnitude  to  prevent  r(nT)  from  overflowing.  Thus,  again 
the  scaling  is  accomplished  by  effectively  reducing  the  filter  gain. 

This  gain  may  be  reintroduced  at  other  points  in  the  system  If  desirable, 
j The  error  analysis  for  the  scaled  filters  or  control  algorithms  follows 

from  previous  results. 

Another  Important  alternative  is  the  use  of  an  adaptive  scale  factor 
that  is  utilized  only  when  an  overflow  is  detected  in  the  system.  If 
this  occurs  and  if  sufficient  computational  time  is  available  the  filter 
can  be  automatically  scaled  until  overflow  ceases  to  exist. 


I 


>1 


i 

•I 


EXAMPLE 


1 


To  Illustrate  the  procedure  consider  the  system  shown  In  Fig.  1 when 

C(z)  is  a discrete  Integral  plus  proportional  controller 

(k,  T/2  + k \z  + k,  T/2  - k 
\ 1 p/  1 £ 

z - 1 


D(z)  = 


(17) 


and  the  plant  Is  given  by  G(s)  = 2/s+5.  The  sampling  Interval  was 
chosen  as  .01  and  the  constants  were  first  assigned  values  = .75 
and  ~ .25.  This  results  in  a sluggish  system  with  no  overshoot  for 
a step  input.  For  a step  input  the  maximum  value  of  e(nT)  is  one 
and  the  maximum  value  of  the  controller  output  r(nT)  is  5/2.  To 
obtain  a controller  with  a gain  of  one,  a multiplying  factor  of  5/2 
is  placed  in  the  data  reconstruction  element.  The  resulting  z transforms 
are 


G’(z)  = 


and 


D’(z)  = 


1 _L. 

2 s s+5 


.1015  z - .0985 


.04877 
Z-. 95123 


(18) 


z - 1 


(19) 


Plots  of  e(nT),  r(nT)  and  y(t)  are  shown  for  the  scaled  system 
with  x(t)  = y(t)  in  Fig.  5. 

For  a direct  form  one  implementation  shown  in  Fig.  6 there,  are 
two  input  rounding  errors.  The  two  product  rounding  errors  are  introduced 
between  the  numerator  and  denominator  portions  of  D'(z).  If  the  errors  are 
each  assumed  to  be  step  functions  with  magnitude  q/2,  the  method  developed 
by  Slaughter  [9]]  results  in  a steady  state  error  at  the  system  output 
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G'(2) 


E = lim  e[nT]  = 
ss 

n-»-» 


/ q 2 G' (z)  D'  (z)  ^ ^ z-1 

\\l  + G'(z)  D'(z)J  [l  + G’(z)  D' 


(z) J J z=l  ’ 


which  gives  upon  substitution 


E = (33A.33)q 
ss 


For  a direct  form  two  implementation  there  are  two  input  rounding 
errors.  The  two  product  roundings  are  Introduced  at  the  output  of 
D'(z).  The  steady  state  error  is  given  by 

P _ r Ti  - / q z G'  (z)  D'  (z)  q z G'  (z)  1 

Ess  ^ ■ I 1 + G'(z)  D'(z)  1 + G'(z)  D'(z)  [ z=l 

n->^  i J 

which  yields 


E = q 
ss 


Substituting  (14)  into  (21)  and  (23)  now  gives  an  expression  for 

the  word  length  in  terms  of  E /A  as 

ss 


^dfl  - ^°®2 


J(334.: 
(E  > 

^ SS 


for  the  direct  form  one  implementation  and 


^df2  - {(E  /A) 

ss 


for  the  direct  form  two  implementation.  Note  that  the  variable  A refers 
to  the  maximum  value  of  the  variable  e(nT)  for  a controller  with  a gain 
of  one  and  no  overflow.  Values  of  N are  tabulated  in  Table  1. 
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1 


If  the  rounding  errors  are  modeled  as  Independent  random  variables 

2 

with  mean  zero  and  variance  q /12,  the  steady  state  variance  of  the 
error  at  the  system  output  for  a direct  form  one  Implementation  Is  given  by 


Var (e) 


12 


{! 

U=o 


h (kX)  + I 
k 


h^(kT)l 

0 J 


(26) 


where  h(kX)  Is  the  Inverse  z transform  of 


H'(z)  = 


D’(z)  G’(z) 

1 + D’ (z)  G' (z) 


(27) 


and  hj^(kX)  Is  the  Inverse  z transfer 


Hi  (z)  = 


G’(z)/(z-l) 

1 + G’(z)  D'(z) 


Substituting  gives 


(28) 


Var(e)  = (25.3)q 


(29) 


For  the  direct  form  two  Implementation 


Var(e)  = (4.07  x 10  ^)q^ 


(30) 


To  establish  word  length  requirements  the  variable  e/A  Is  now 
assumed  to  be  a Gaussian  random  variable.  The  probability  that  e/A  Is 
between  where  Is  the  previously  given  steady  state  error. 


Is  now  defined  as 


Prob 


{- 


£ E 

SS  £ ss 

A.  - A - A 


confidence  level. 


(31) 
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For  a given  confidence  level,  normalized  distribution  tables  can  be 


used  by  writing 


!ss  ^ ^ ^ fss  1 
1 Ao  — Ao  — Ao  J 


confidence  level 


where 


2^-1 


direct  form  one 


^•07  X lO"^ 


2^-1 


direct  form  two 


For  a 95%  confidence  level 


and  substitution  of  (32)  into  (33)  now  gives 


^dfl  - ^°^2  If.  /a 


for  the  direct  form  one.  For  the  direct  form  two 


^df2  - ^°^2  1 E /A  ^ 


The  value  of  N,.,  and  N,,.„  are  also  tabulated  in  Table  1. 
at  i at  z 

Referring  to  Table  1,  the  form  that  is  used  to  Implement  the  controller 
is  obviously  a critical  factor.  The  use  of  the  variance  of  error  results 


in  a somewhat  shorter  word  length  for  a given  value  of  E /A.  However, 

ss 

in  practice  a very  short  word  length  will  give  rise  to  other  undesirable 
effects  due  to  large  deadbands  caused  by  the  rounding.  Therefore,  for  the 
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example,  a word  length  of  ten  bits  plus  sign  might  be  utilized.  This 

_3 

gives  an  error  E^^/A  less  than  10  for  the  controller  implemented 
in  direct  form  two. 

To  improve  the  response  a new  control  algorithm  was  designed 

using  K = 3 and  K.  = 15.  The  resulting  discrete  controller  transfer 
P i 

function  is 


D’tz]  = ~ . 

z - 1 

Simulation,  however,  shows  the  controller  output  will  overflow  so  a value 
of  1.23  was  factored  from  (36)  giving 


(36) 


D"[2] 


z - .9593495935 
z - 1 


and 


G"[z] 


.0599878078 
z - .9512294245 


(37) 


(38) 


Response  curves  for  the  resulting  system  are  shown  in  Fig.  6.  Note 
that  the  steady  state  output  for  r(nT)  is  about  .8  units  which 
means  the  full  range  of  the  digital  controller  is  not  being  utilized 
fully.  Word  length  requirements  for  the  system  are  given  in  Table  2. 

For  the  direct  form  one  the  improved  system  requires  fewer  bits.  While 
this  result  may  not  be  expected  careful  examination  of  the  system 
structure  shows  the  new  system  to  be  less  critical  to  noise  generated 
within  the  system  when  the  controller  is  implemented  In  direct  form  one. 
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TABLE  1 


TABLE  OF  WORD  LENGTH  REQUIREMENTS 

= .75  K = .25 
1 P 
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TABLE  2 

TABLE  OF  WORD  LENGTH  REQUIREMENTS 

K = 15  K = 3 
i P 
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APPENDIX  A 


In  this  section  a closed  form  method  for  computing  the  discrete 
autocovariance  function  is  derived.  Consider  a discrete  transfer 
function 


H[z]  = 


M 

I a z 
k=0 


-k 


N 


1 - y b.  z“^ 


= Hj^[z]  H^fz] 


j = l ^ 


where 


H^[z]  = 


1 - y b , z ^ 


j=i 


and 


H,[z]  = y 


2 k=0 


-k 


If  an  independent  noise  sequence  {e[jT]}  with  mean  zero  and  variance 


o is  applied  to  H^^[z]  the  difference  equation  representation  is 


given  by 


N 


x[nT]  = lb  x[(n-j)T]  + e[nT] 

j=l  ^ 


The  discrete  autovariance  at  time  t = nT  is  now  defined  as 


(la) 


(2a) 


(3a) 


(4a) 


R [kT]  = E{x[nT]  x[(n-k)T]}  k = 0,  1,  2,  . . .n  . (5a) 

n 


Since  e[nT]  is  taken  from  an  independent  noise  sequence  x[(n-k)T] 
and  e[nT]  are  uncorrelated.  Multiplying  (4a)  by  x[(n-k)T], 
k = 0,  1,  2,  ...  N-1  and  taking  the  expectation  now  results  in  a set  of 
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coupled  difference  equations  in  the  variables  {R  [OT] , R [T]  R [(N-l)T]} 

n n n 

Taking  the  2 transform  and  applying  the  final  value  theorm  now  permits 
the  first  N terms  of  the  steady  state  autovariance  sequence  to  be 
written  in  matrix  form  as 


""  R[0T] 

r 2'i 

a 

R[T] 

0 

. 

= 

0 

0 

_R[(n-l)T]_ 

0 _ 

where  M is  the  matrix 


N a 

N-1 

N-2 

V ,2 

^ 1 
3=1  ^ 

(b^ 

(1  + b2) 

• 

• 

• 

(b,) 

^^N-2 

• 

1 


J 


The  remaining  terms  in  the  autovariance  sequence  are  found  from 

N 

R[kT]  = - [ b.  R[(k-j)T]  k N . 

j=l  ^ 

The  steady  state  autovariance  at  the  output  of  H[z]  is  given  by 
11m  E{y[nT]  y[(n-k)T]}  = R [kT]  = a 


c a 


c = R[(k-j+l)T]  = 11m  E[x[nT]  x[n-(k-j+i)T]  . (10a) 

J n-^ 


When  computing  the  steady  state  autovariance  for  k < M-1  certain  terms 
In  the  C matrix  are  found  using  the  relationship  R[kT]  = R[-kT). 

Thus,  a fixed  nonrecursive  expression  only  exists  for  k ^ M-1. 

The  variance  at  the  output  of  the  system  is  found  using  the  results 
of  (6a) . The  expression  at  the  output  is 

Tim  E[v  = a*"  C,  a (11a) 

■ n 1 

where  is  the  autovariance  matrix 
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in  Direct  Form 


CVJ 


a. 


TD 

C 

o 


o 

II 


</) 


c 


U) 

0) 

> 


3 

o 


^n 

I 


0) 

or 


40 

0) 


3 

9* 

u. 


apniiidiuv 


CJ 

o 
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ABSTRACT 


A strategy  for  quantizing  the  coefficients  in  a general  second  order 
digital  control  algorithm  is  presented.  First,  the  errors  in  the  magnitude 
and  phase  functions  for  the  algorithm  are  derived  in  terms  of  the  filter 
coefficients.  By  specifying  a maximum  allowable  error  for  each  function 
over  a given  frequency  range , quantization  regions  can  then  be  established. 
An  example  consisting  of  a lag-lead  digital  control  algorithm  is  Included 
to  Illustrate  the  procedure. 


INTRODUCTION 


The  use  of  discrete  data  processing  algorithms  as  an  integral  part 
of  a control  system  is  well  known.  Books  to  Tou  [1],  Monroe  [2]  and  Kuo  [3] 
describe  many  of  the  design  concepts  for  these  systems  called  sampled 
data  or  digital  control  systems.  In  a somewhat  independent  effort,  design 
and  implementation  procedures  have  been  developed  for  digital  filters. 

A collection  of  papers  in  this  area  has  recently  been  edited  by  Rabiner 
and  Rader  [4]. 

The  accuracy  of  any  digital  control  algorithm  is  dependent  on  the 
finite  word  length  used  to  represent  the  input  data,  the  filter  coefficients 
and  the  intermediate  products.  A recent  review  paper  on  this  topic  that 
is  concerned  with  open  loop  digital  filters  is  by  Liu  [5]. 

Error  due  to  deterministic  coefficients  in  digital  filters  is  described 
by  Kaiser  [6],  Otnes  and  McNamee  [7]  and  Mantey  [8].  Gold  and  Rader  [9,10] 
have  investigated  the  effect  of  coefficient  error  on  the  filter  pole 
locations  and  have  developed  filter  forms  that  are  less  sensitive  to  error. 

Avenhaus  [11]  developed  a procedure  for  finding  the  coefficient  word  length 
in  a digital  filter  when  rounding  is  employed.  This  procedure  assumes 
an  ideal  magnitude  function  in  the  error  formulation  and  results  are 
developed  only  for  the  passband  and  stopband.  A statistical  approach  to 

the  problem  is  described  by  Knowles  and  Olcayto  [12].  7 

1 

A general  result  from  the  references  is  that  higher  order  digital 
filters  and/or  control  algorithms  should  be  implemented  as  a cascade  or  j 

parallel  connection  of  first  and  second  order  filter  sections.  Wlien  the  | 

filter  sections  have  been  determined,  the  filter  coefficients  must  then  j 

be  quantized  in  such  a way  that  the  specifications  are  satisfied.  1 
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Tile  implementat ion  of  digital  control  algorithms  with  current  digital 
hardware  can  also  place  constraints  on  the  coefficient  quantization.  Using 
medium  or  large  scale  Integrated  circuits  there  are  often  constraints  on 
the  input  and  output  pin  connections.  Computations  done  within  the 
integrated  circuit  on  the  other  hand  can  often  utilize  more  bits.  This 
means,  for  example,  that  sums  of  products  can  often  be  truncated  only 
after  all  multiply-add  operations  have  been  completed.  Since  the 
coefficients  represent  a set  of  inputs  to  the  control  algorithm  quantiza- 
tion should  be  carried  out  that  allows  a minimum  number  of  bits  to 
represent  each  coefficient. 

When  the  control  algorithm  is  implemented  using  fixed  point  arithmetic 
on  a computer  with  a word  length  of  n bits,  it  is  desirable  to  quantize 
the  coefficients  to  less  than  n bits.  This  allows  flexibility  in 
scaling  when  the  programming  is  carried  out. 

In  this  paper  a strategy  for  quantizing  the  coefficients  in  a second 
order  control  algorithm  is  developed.  First,  the  errors  in  the  magnitude 
and  phase  functions  are  expressed  in  terms  of  the  filter  coefficients. 

By  comparing  the  derivatives  of  the  magnitude  and  phase  functions  with 
respect  to  each  coefficient,  closed  regions  can  be  established  for  the 
quantization  of  the  coefficients  over  a frequency  range  of  interest.  These 
regions  are  also  useful  when  the  coefficients  are  quantized  using  conventional 
rounding.  The  procedure  is  illustrated  for  a second  order  lag-lead 
compensator  designed  using  the  bilinear  z transform. 
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ERROR  ANALYSIS 


Magnitude  Function 


Given  a second  order  discrete  transfer  function 


H(z)  = 


^0  " + z + a^ 

z^  + z + b2 


the  frequency  function  is  given  by  letting  z = e^^^.  The  resulting 
magnitude  function  is  written  as 

l“<“>i  ■ 

where 


(2  2 2 

A(w)  = ^ ^2  ^ ^^1^2^  coswT  + 23^3^  cos2wT 


B(w)  = /l  + + 2b^(l  + b2)cosuT  + 2b2CosuT|' 


When  the  coefficients  3re  allowed  to  vary  (2)  can  be  considered  as  a 
function  of  five  variables  for  any  value  of  frequency.  This  leads  to 
the  differential  approximation 


a1h(w)|  = 


Ab^  + 


Aao  + 


Aa^^  + 


where  each  partial  derivative  is  a function  of  frequency  given  by 


- A(o)){bj^  + (1  + b2)cosajT} 


i 

J 
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(7) 


and 


siH(^) ! ^ 
3b^ 


A(u){b2  + b^  cosojT  + cos2u)T} 


3|H(ai)  I ^ 

3a„ 


cosojT  + cos2wT 
A(ui)  B(o)) 


(8) 


3|H(oo)|  _ 


3a, 


A(w)  B(w) 


(9) 


3 |H(ui) 
3a_ 


a.  + a,  coswT  + a_  cos2wT 


A(o))  B(u) 


(10) 


In  the  above  expressions  A(w)  and  B(w)  are  given  by  (3)  and  (4) 


Phase  Function 

Given  the  transfer  function  described  by  (1)  the  resulting  phase 
function  is 


6(u)  = Tan 


j a^  sln2u)T  + sinooT 
a^  cos2u)T  + a^  cosuT  + a^ 


Tan 


-1 


r sln2u)T  + b^  sini  jT  'I 

lcos2(i>T  + bj^  cosuT  + b2  J 


For  simplicity  let 


Nj^(io)  = 3q  cos2wT  + sintoT  + a^ 

N2(w)  = a^  sln2wT  + sinuT 

Dj^(uj)  = cos2uT  + bj^  coaoiT  + b2 

and 

D2(uj)  = sln2wT  + b^  slnoaT 


(11) 


(12) 

(13) 

(14) 

(15) 
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Th'i  differential  phase  approximation  is  now  described  by 
^ 36(w)  . 36(u)  , 39((;j)  . 


where 


. 50  , 9e  , 

+ r — Aa  + v;—  Aa„ 

3a^  1 da^  2 


39(gj) 


-|di(cj)  sinojT  - coswt| 


3b^ 

(00)]^  + [D2(g))]^ 

39(g)) 

02(1^) 

3b2 

[D^(g)) 

+ [0,(0))]^  ’ 

39(y) 

N^(g)) 

sin2g)T  - N2(g))  cos2g)T 

[N^( 

g))]^  + [N2(g))]^ 

39(g)) 

N^(g)) 

sing)T  - N2(g))  sing)T 

[N^(co)]^  + [N^Cw)]^ 


39(g)) 


- N2(gj) 


aa^  [N^(gj)]^  + [N^Cg))]^ 


Determination  of  Quantization  Regions 

The  expressions  described  by  (5)  and  (16)  are  linear  algebraic 
equations  in  the  variables  Ab^,  Ab^,  Aa^,  Aa^^  and  Aa2  with  coefficients 
that  are  functions  of  frequency.  For  specifications  described  by  the 
bounds  “Ej^  j1  a|  11(g))  | £ and  2 — A9(g))  ^ E2>  (5)  and  (16)  each  form 

a five  dimensional  space  and  the  strategy  is  to  first  find  the  two  smallest 
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spaces  by  utilizing  the  proper  combination  of  coefficients.  The  vector 


(Ab^  Ab^  A3q  Aa^  Aa^ } is  then  forced  to  be  inside,  of  each  space 
by  the  proper  quantization  of  eacli  filter  coefficient. 

In  practice,  however,  the  quantization  problem  can  often  be  treated 
using  two  dimensional  regions.  In  the  example  that  follows,  one  two 
dimensional  region  is  sufficient  for  the  specification  on  AjH(ui)]  and 
one  region  for  the  specification  on  A6(aj).  In  other  practical  algorithms, 
the  error  in  AH(aj)  and  A6(u)  can  be  distributed  among  the  coefficients 
so  that  one  and  two  dimension  regions  can  be  utilized  for  the  actual 
quantization. 


E.XAMPLE 


For  an  analog  lag-lead  compensator  described  by 


C(s) 


(s  + . 1)  (s  + 1) 
(s  + .01)  (s  + 10) 


the  equivalent  discrete  compensator  designed  using  tiie  extended 
z transform  is  given  by 

. . 0.8356618816z^  - 1.626383584z  + 0.7909250553 

c(z)  = . 

z - 1.592691562Z  + 0.592894916 


bilinear 


(22) 


(23) 


In  the  design,  the  sampling  time  was  chosen  to  be  T = 0.05. 

Graphs  of  tlic  partial  derivatives  described  by  (6)  througli  (10)  and 
(17)  through  (21)  are  shown  in  Figs.  1 through  5 as  functions  of  frequency. 
Values  for  eacli  partial  derivative  were  computed  over  a range  of  zero  to 
twenty  radians,  however,  to  illustrate  the  important  features  the  frequency 
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range  for  each  derivative  was  reduced.  For  higher  frequencies  each  partial 
derivative  asymptotically  approaches  zero. 

Referring  to  the  graphs  the  following  assumptions  can  be  made.  First, 
equations  (6)  and  (7)  are  approximately  the  same  as  are  (17)  and  (18). 
Secondly,  equations  (8),  (9)  and  (10)  are  approximately  the  same  as  are 
(19),  (20)  and  (21).  From  computer  printout,  the  approximations  are  good 
to  three  decimal  places.  This  Implies  that  (5)  and  (16)  can  be  rewritten  as 


A|H(o)) 


+ Aa^  } 


+ 


+ Ab^} 


(24) 


Ae(u)  = ^ {Aa_  + Aa  + Aa  } 

da^  0 12 


+ {Ab  + Ab  } 

3bi 


Quantization  regions  are  now  established  by  letting  Ax^  = {Aa^  + Aa^  + Aa^l 
and  Ax^  = (Ab^  + Ab2)  and  specifying  values  of  A|H(aj)j  and  A0(w).  For 
A|h(uj)|  = *0.1  the  minimum  quantization  region  using  (24)  is  shown  in 
Fig.  6.  This  region  is  established  by  substituting  various  values  of 
0|H(oj)|/9aQ  and  3 |h(uj)  | / 3b  into  (24)  to  find  the  set  of  lines  that 
enclose  a minimum  area.  For  0(a))  = *1  degree  the  minimum  region  is 
shown  in  Fig.  7.  The  actual  coefficient  quantization  is  done  so  that  Ax^^ 
and  AX2  lie  Inside  of  each  quantization  region.  The  quantized  coefficients 
and  quantization  errors  are  shown  in  Table  1.  Note  that  the  quantization 


i 


148 


i 


was  done  to  yield  coefficients  that  can  be  represented  by  a minimum  number 
of  binary  bits.  The  binary  two's  complement  representation  for  each 
coefficient  is  also  given  in  Table  1.  Graphs  of  the  error  plotted  as  a 
function  of  frequency  are  shown  in  Figs.  8 and  9.  These  graphs  show  that 
the  specifications  are  met.  For  higher  frequencies  the  errors  converge 
toward  zero. 

The  reader  should  be  aware  that  quantization  leading  to  a binary 
representation  for  each  coefficient  is  more  difficult  than  quantizing 
a decimal  number.  To  achieve  a minimum  word  length  it  is  necessary  to 
carefully  quantize  each  group  of  coefficients  so  that  error  cancellation 
takes  place.  The  results  obtained  in  Table  1 are  now  compared  to  those 
obtained  using  conventional  rounding.  First  the  coefficients  in  the  lag- 
lead  compensator  were  converted  to  18  bits  plus  sign  as  shown  in  Table  2.  | 

Then  conventional  rounding  was  used  to  obtain  fifteen  bit  plus  sign  to 
12  bit  plus  sign  representation  for  the  coefficients.  Rounding  was  done 
by  looking  at  the  most  significant  bit  of  the  truncated  portion  of  the 
coefficient.  If  this  bit  is  a one,  a one  is  added  to  the  least  significant 
bit  of  the  truncated  coefficient.  If  the  bit  is  zero  a zero  is  added  to 
the  least  significant  bit.  The  truncated  coefficients  are  shown  in 
Table  3.  The  magnitude  function  error  curves  are  shown  in  Fig.  10.  Here, 
with  conventional  rounding  at  least  lA  bits  plus  sign  are  needed  to 

1 

satisfy  the  specification  a|H(o))|  = + .1.  Furthermore,  for  this  ^ 

example  when  12  bits  plus  sign  are  used,  conventional  rounding  gives 

coefficient  values  tliat  result  in  a dc  gain  of  0.48828042  x 10^  rather  j 

than  the  desired  dc  gain  of  one.  Wlien  the  proper  strategy  is  used,  liowever, 

onlv  12  bits  plus  sign  are  needed  for  the  coefficients  as  shown  in  Table  1.  j 
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10.01101000010001011 
b,  0.100101111100011110 

0.110101011110110110 
10.01011111101001011 
a^  0.110010100111101000 

TABLE  2 

EIGHTEEN  BIT  PLUS  SIGN  BINARY  EQUIVALENTS 
FOR  THE  FILTER  GOEFFICIENTS 
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10.01101000010001 

Abi 

= 2.08  X 10“^ 

*^2 

0.100101111100100 

Ab2 

= -5.92  X lO”^ 

^0 

0.110101011110111 

Aao 

= -9.60  X lO"^ 

^1 

10.01011111101001 

Aa^ 

= 2.02  X 10*^ 
-7 

^2 

0.110010100111101 

Aa^ 

= 9.83  X 10 

(a) 

-5 

AX2  = 2.02  X 10 

Fifteen  Bits  Plus  Sign 

Ax^  = 2.02  X 10 

•’l 

10.0110100001001 

Ab  = -4.02  X 10~ 

'■  -7 

'^2 

0.10010111110010 

Ab  = -5.92  X 10 
^ -5 

^0 

0.11010101111011 

Aa  = 2.96  X 10 
^ -5 

10.0101111110101 

Aa  = -4.08  X 10 
^ ^-5 

^2 

0.11001010011111 

Aa^j  = -2.95  X 10 

(b) 

-5 

Fourteen  Bits  Plus  Sign 

AXj^  = -4.07  X 10  AX2  = -4.08 

10.011010000100 

Ab^  = 

8.19  X 10' 

^-7 

0.1001011111001 

Ab2  = 

-5.92  X 10 

„-5 

®0 

0.1101010111110 

II 

0 

< 

-3.15  X 10 
-5 

^1 

10.010111111010 

Aa^  = 

8.13  X 10 

-5 

®2 

0. 1100101001111 

Aa2  = 

3.15  X 10 

(c) 

Thirteen  Bits  Plus  Sign 

Ax^  = 8.13  X 10”^ 

Ax^  = 8.13  X 

‘’I 

10.01101000010 

Abi 

= 8.19  X 10” 

„-5 

0.100101111100 

Ab2 

= 12.15  X 10 

-5 

^0 

0.110101011111 

Aao 

= -3.15  X 10 
-5 

“1 

10.01011111101 

Aa^ 

= 8.13  X 10 

„-5 

^2 

0. 110010101000 

Aa2 

=■  -9.06  X 10 

(d) 

Twelve  Bits  Plus  Sign 

AXj^  = -4.08  X 10 

AX2  = 20.34 

Table  3 

Filter  Coefficients  tibtained  Usinn 
Conventional  Rounding 
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